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SUMMARY 


A steady-state design computer program has been developed to predict the 
performance of pumping rings as functions of geometry, applied loading, speed, 
ring modulus, and fluid viscosity. Additional analyses have been developed to 
predict transient behavior of the ring and the effects of temperature rises 
occurring in the hydrodynamic film between the ring and shaft. The analysis 
was initially compared with previous experimental data and then used to design 
additional rings for further testing. 

Tests were performed with Rulon, carbon-graphite, and babbitt rings. Two 
different shaft diameters were used for the babbitt rings. The design analy- 
sis was used to size all of the rings and to select the ranges of clearances, 
thickness, and loading. Although full quantitative agreement was lacking, 
relative agreement existed in that rings that were predicted to perform well 
theoretically, generally performed well experimentally. Some causes for 
discrepancies between theory and experiment are believed to be due to starva- 
tion, leakage past the secondary seal at high pressures, and uncertainties in 
the small clearances and local inlet temperatures to the pumping ring. 

The design criteria that evolved require the applied loading to be of the 
order of the desired pumping pressure, the flow requirements and tolerances to 
dictate clearance, and the elastic modulus and ring compliance to be such that 
the deflection under load statically results in clamping at very small inter- 
ferences so that back flow is inhibited, but that excessive power loss and 
wear do not occur. 

It was found that the pumping ring could be used to generate its own loading 
pressure without any priming if an initial taper was applied to the ring. 
However, for untapered rings, an initial loading had to be applied before 
self-pumping could be obtained. 

A separate preliminary analysis has been performed for a pumping Leningrader 
seal. This analysis can be used to predict the film thickness and flow rate 
through the seal as a function of pressure, speed, loading, and geometry. 
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1.0 INTRODUCTION 


An analysis of pumping rings was performed under certain simplifying condi- 
tions in previous phases of this work [ 1 ] . These conditions consisted of 
first ignoring the contribution, if any, of the back flow occurring during the 
return stroke of the rod.’ Other simplifications related to the use of 
constant or average parameters, namely constant viscosity and mean rod veloci- 
ty, throughout the stroke. The latter also presumed the neglect of squeeze 
film forces due to the variation of film thickness engendered by the harmonic 
motion of the rod. 

The comparison of experimental data with theoretical results based on the 
simplified analysis showed good agreement with respect to maximum pressures 
generated by the pumping ring. The flows produced at reservoir pressures 
below the maximum, however, were consistently lower in the experiments than 
those indicated by the analysis. The agreement between theory and experiment 
for the carbon-graphite rings [l] has been found to be incorrect due to an 
erroneous use of an excessively low viscosity in the theoretical computation. 

The present work, an extension of the previous effort, is aimed at advancing 
the analysis of pumping rings. Thermal effects and variable rod velocity were 
included and, with it, the effects of squeeze film action in the fluid film. 
The analysis was extended to include the backstroke and concurrent cavitation 
and their effect on net flow. The new analysis was then used to run a parame- 
tric study in order to obtain optimized configurations of pumping rings of 
different shapes and materials including the effects of' a geometric taper. 
The results of tests run on these optimized designs were then compared with 
calculations based on the new analysis. 

A separate, preliminary anaLysis of the pumping Leningrader seal has also been 
performed. . Since this analysis is separate from the pumping ring work, its 
results are presented in Appendix A. 
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2.0 ANALYSIS WITH CONSTANT PARAMETERS 


2.1 Basic Equations 

The equation governing deflection, w, of an axisymmetric shell under bending 
is . 


12(l-v 2 ) 


4 

d w Et . 

— r + — r w = -P(x) 

dx R 


(2-1) 


Referring to Figure 2-1, the radially outward loading, P(x), may be expressed in 
terms of the clamping load, p Q , and the hydrodynamic pressure, p, as follows: 

{ 0, -L x £ x £ 0 

p, 0<x£L-e (2-2) 

p-p Q , L - e < x 1 L 


The hydrodynamic pressure, p, is determined from the solution of the Reynolds 
Equation 


= 

dx 


6 yU 

o 



(2-3) 


with the geometry of the hydrodynamic film as given in Figure 2-2. k is a 
constant of integration related to the flow, and U q is the average speed, which 
is related to the frequency, f, and the stroke, s, by: 


U = 2fs 
o 

’ » 

This velocity represents the average of the sinusoidal velocity over each half 
cycle. If the system contains two opposing pumping rings, it is double acting 
and the average velocity, U q is assumed to prevail over the entire cycle (though 
in two different pumping rings) for the forward stroke, as well as for the back- 
stroke. In this manner, the transient problem is reduced to and simplified into 
a quasi-steady-state process. For a single ring, the resulting flow should be 
divided by 2. 
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The film thickness, h, appearing in Equation (2-3) is related to the deflection, 
w, by 

h=(C-w),0<x<L 


The system of equations given by Equations (2-1) through (2-3) represents a 
fifth-order set of differential equations requiring five boundary conditions, 
in addition to a sixth condition for the evaluation of the constant, k, appear- 
ing in Equation (2-3) . Two conditions result from the clamped-end requirement 


dw . 

w =. -r- = 0 at 
dx 


x = -L, 


Two conditions resulting from the prescribed pressures at x = 0 and L 


p = 0 at x = 0 p = p^ at x = 1 


where p^ is the sealed pressure. The remaining two conditions relate to the 
free-end requiring zero moment and zero shear, or 




0 at x = L 


The method of solving this set of equations subject to the specified boundary 
conditions is outlined in Reference [ 1] . The solution and the results of this 
analysis, reported in Reference [1], are based on the constancy of both viscosi- 
ty and speed, namely 

y = y = constant 
o 

U = U = 2fs = constant 
o 


The analysis in Reference [1] did not consider the backstroke and was thus only 
applicable to sufficiently high loads and low elastic moduli to result in clamp- 
ing during the backstroke. 
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2.2 Simplified Approach 


A scrutiny of the pumping ring solutions formulated in Section 2.1 shows that 
the shape of the film over the hydrodynamic portion of the ring is nearly 
linear for the entire range of ring parameters and operating conditions. A 
few such examples are shown in Figure 2-3. Consequently, the problem can be 
considerably simplified if one postulates that the configuration of the film 
is tapered similar to that of a plane slider. In addition, a constant taper 
(shown in Figure 2-2) makes it possible to later treat the backstroke and 
accompanying cavitation. Mathematically a constant taper implies that h' = 
constant. The pertinent expressions for the hydrodynamic component of pumping 
ring action are then given by 

h (O = h 2 + Ah (1 - O (2-4) 

p(h) = l/(h 2 - h x ) - 1/h + k'/(2h 2 ) + C ± (2-5) 

where 


k' = [2h 1 h 2 /(h 1 +h 2 )j[l - (h 1 h 2 p f )] 


c x » [i/(h 1 h 2 ) ] [1 + (h 2 p f )] 


The elastic deformation equation remains the same and may be written in dimen- 
sionless form as 

Cl for -L < l < 0 

o (d 4 h/d? 4 ) + h= jl + BpforO<£<l-e • (2-6) 

L 1 + 3(p - p Q ) for 1 - t < £ < 1 

The solution algorithm thus consists of determining values for h 2 and Ah by 
the use of the secant method. Equation (2-6) is solved subject to the 
previously stated elasticity boundary conditions. Convergence is achieved 
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when the values of h(l) and h'(l) computed from the solution to Equation (2-6) 
are within the prescribed tolerance limit of the values of h 2 and -Afi used in 
calculating p from Equations (2-4) and (2-5). 

2.3 The Backstroke 


2.3.1 Analytical Approach 

The previous analysis assumed clamping during the backstroke. This is valid 
for very high loadings, or for highly flexible pumping rings. Lower ring 
loadings which do not cause clamping during the entire backstroke may be 
desirable for pumping ring design in order to reduce wear. Also, when the 
upstream pressure is high, the ring may stay open during the reverse stroke, 
even under high clamping forces. Thus, backstroke effects are here added on 
to the analytical model. 

The basic equations remain unchanged except for two important aspects. One is 
that the shaft motion will be in the negative x direction; thus Reynolds 
Equation assumes the form 

dp/d? = -[(h - K)/h 3 ] (2-7) 


The other critical modification consists in the appearance of cavitation. As 
shown in Figure 2-4, due to the divergence of the film in the direction of 
motion, there may not be enough fluid to fill the gap at sufficiently high 
values of h. The fluid film will then break up into a pattern of streamers 
similar to that which occurs in the diverging films of a journal bearing. 
From continuity requirements, the downstream boundary conditions of ? c , where 
the film ends, must then satisfy the following boundary conditions: 

(dp/d?) at ? = £ =0 

c 

( 2 - 8 ) 

p(S c ) = 0 
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To account for cavitation during the backstroke, the equations are integrated 
backward from £ = 1 . As for the forward stroke solution, values of h2 and Ah 
must be determined by using the secant method. Selected values for h2 and Ah 
are used in determining h(5) from Equation (2-4) which is in turn substituted 
in the Reynolds Equation (2-7). The resulting equation is then integrated 
analytically, subject to the constraints that p(l) = , and either the condi- 

tion given by Equation (2-8) for 0 < 5 C < 1 or the condition p(0) =;0 if cavi- 
tation is not predicted to occur. These conditions are sufficient to 
determine the constant, K, C c (if applicable), and the pressure distribution, 
p(5). Equation (2-6) may now be integrated backward from 5 = 1» with the 
boundary condition at £ = 1 such that 

h( 1 ) = h2 

h’(l) = -Ah 

h"(l) = h" ' (1) = 0 

The secant method is then used to find the values of h2 and Ah that make the 
quantities |h(-Li) - 1| and |h'(-Li)| be within prescribed tolerance limits. 
The solutions so obtained provide dimensionless values of the film thickness 
profile, the pressure profiles and the flow rates for prescribed values of 
parameters Ot, 8, p q , p^ , £, and Lj. 

2.3.2 Nature of Solution 

The nature of the solution for the backstroke depends very much on the value 
of p£ relative to the clamping force, p Q . It is thus first necessary to 
describe the conditions which are apt to generate either high or low levels of 
Pf . 

If a pumping ring delivers oil to a closed reservoir of finite volume, the 

parameters a, 8^ E, and p describing the pumping ring will predetermine the 

o 

maximum p^ that the pumping ring is likely to generate. At that point, i.e., 

when p is reached, there will be no further inflow into the reservoir, and 
fm 

the pressure in the reservoir will be maintained by pumping ring action at 

p . In Figure 2-5 this situation would be represented by both valves A and B 
fm 

being closed. However, should there be an outflow from the reservoir, that 
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Fig. 2-5 Conditions Determining Level of 
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is, were vaLve A open, then the pumping ring would be delivering a net inflow 
equal to the amount of outflow, with the pressure in the reservoir below the 
maximum, i.e., p^ < . On the other hand, if an external pump were to deliv- 

er a certain flow, then it is possible to have p^ > Pf m * The pumping ring 
would then be overloaded. The net flow over a cycle would be negative and the 
outflow through the pumping ring would equal the inflow delivered by the pump. 
Since there is no external pump in the present system, the overloaded condi- 
tion will be of academic interest only. 

The interaction of pumping ring behavior versus conditions in the reservoir 
has a direct bearing on the nature of the backstroke solution. The mechanism 
can be summarized as follows: 

• If Pf is low, i.e., if there is an appreciable outflow from the reser- 
voir, the film during the backstroke always cavitates. It is only by 
having a cavitating backstroke that a net inflow can be maintained over 
a complete cycle. 

• At high values of p^ , there are generally two possible solutions, name- 
ly: 

- High h2 with a noncavitating film 

- Low h2 with a cavitating film 

V:' ■ 

Such a double set of possible solutions is shown in Figure 2-6, both 
the same pj = 2.7 and p Q = 3.5. The performance of the pumping ring 
under these two sets of conditions is as follows: 



Noncavitating 

Cavitat ing 

Back h2 

0.482 

0.120 

Forward h2 

, 0.580 

0.580 

Back K 

-1.433 

-0.225 

Forward K 

-0.42 

-0.42 

Net K 

-1.85 

-0.642 

^cav 

— 

0.87 
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Solutions for High Values of pf 





• If there is no supply of oil to the reservoir by a pump, that is, if the 
pressure, p^, is generated by the pumping ring only, then the solution 
of a high h2 with a noncavitating film is practically impossible, since 
a thick noncavitating film would produce a net outflow so that no high 
p^ producing that outflow could possibly be generated in the first 
place. 

Thus, in conclusion, for practical cases of an active pumping ring, the back- 
stroke is always accompanied by cavitation. 

2.3.3 Performance Characteristics 

For the practical range of ring operation, the pumping ring will, in most 
cases, cavitate during the backstroke. When the condition. of (dp/d£)l£ c = 0 
is introduced during the backstroke, the resulting film thickness, pressures, 
and flows are significantly altered from the full film case. Table 2-1 gives 
the detailed characteristics of cavitating pumping rings as a function of 
sealed pressure, p^. Particularly significant for these purposes is the 
amount of back flow and its effect on the net pumping accomplished over a 
complete cycle. Figures 2-7 and 2-8 show the pressure distribution during the 
forward and reverse strokes for a range of p^ from 0 to 5. As seen, cases of 0 
< p f < 3 produce cavitation, and the fluid film and corresponding positive 
pressures for these cases occupy merely a fraction of the interspace. This 
reduced pressure profile naturaLly has a strong effect on the film thickness 
distribution shown in Figures 2-9 and 2-10. Although h2 for the noncavitating 
cases, pj > 3, is about the same for the forward and backstrokes, there is a 
manyfold (5 to 10 times) decrease in F 12 during a cavitating backstroke. The 
shapes of the pressure profiles for a range of values of p Q (all previous 
•plots were for p Q = 3.0) are shown in Figure 2-11. 

A scrutiny of the data contained in Table 2-1, as well as of the accompanying 
plots, reveals the following interesting features regarding cavitating versus 
noncavitating films. 
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TABLE 2-1 


PERFORMANCE OF CAVITATING PUMPING RINGS 


-•o = 0.0282 8 = 0.257 e = 0.518 p = 3.5 L = 1.57 


Pf 

Flow , K 

^ 

h 

min 


Forward 

Stroke 

Reverse 

Stroke 

Net 

Forward 

Stroke 

Reverse 

Stroke 

0 

0.355 

0.0452 

0.310 

1.0 

0.213 

0.045 

0.74 

1.0 

0.315 

0.0635 

0.252 

0.983 

0.285 

0.045 

1.011 

1.5 

0.244 

0.0759 

0.168 

0.975 

0.342 

0.052 

1.0 

2.0 

0.0877 

0.0954 

-0.008 

0.963 

0.422 

0.059 

1.0 

2.5 

-0.227 

0.903 

-1.130 

0.0604 

0.527 

0.342 

1.0 

3.0 

-0.786 

2.068 

-2.854 

None 

0.652 

0.603 

1.0 

3.5 

-1.6650 

3.284 

-4.949 

None 

0.790 

0.771 

1.0 
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for Reverse Stroke 









• The change from a cavitating to a noncavitating film is rapid. This is 
shown in Figures 2-12 and 2-13 where the transition region 
(2.0 < p£ < 3.0) is plotted in some detail. The onset of cavitation 
seems to occur at about p»£ = 2.5 where £ c - 0. 

• Noncavitating films seem to be accompanied by a net reverse flow even 
during the forward stroke, emphasizing the previous remark that such 
cases would be of little practical interest in the application of pump- 
ing rings. 

• The total net flow, including both the forward and backstrokes, becomes 
negative, even for cavitating films at about p>£ “ 2.0. 
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Fig. 2-12 Film Thickness at Transition from Cavi 
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Fig. 2-13 Pressure Profile at Transition 


XJ * 


t ; 

i- 

• i 

i j 



~ > ■ . I I l ■■■!■!■ i 

0.6 0.7 0.8 0.9 1.0 


rom Cavitating to Noncavitating Condition 





3.0 ANALYSIS WITH VARIABLE PARAMETERS 


This part of the analysis not only abandons the assumptions of constant speed 
and constant viscosity but also considers some other factors which were previ- 
ously neglected. Specifically, the analysis included the following elements: 

• Thermal effects 

• Variable speed 

• Squeeze film action 

• Starvation 

• Nonparallel contours. 

Of course, the analysis also retained the backstroke and accompanying 
cavitation. However, since for the cavitating backstroke, the film thickness is 
of limited extent, the refinements of variable temperature were not included for 
that part of the cycle. Thermal effects and transient effects were considered 
separately to obtain their individual influence. This provides the quantitative 
corrections associated with these effects without adding the high degree of 
complexity of a fully coupled analysis. 

3.1 Thermal Effects 

3.1.1 The Energy Equation 

Since no side leakage exists, the one-dimensional energy equation could be used. 
This, of course, assumes that temperature variation with y can be averaged, a 
fact which, as will be seen later, is particularly relevant here. It was also 
assumed that all the heat generated is convected away by the lubricant. 

Ignoring conduction, the one-dimensional energy equation can be written as, 

pcu(3T/3x) = y[ (3u/9y) 2 ] (3-1) 

Since, from one-dimensional bearing theory 

u - h ^ y <y • h) + u o < H JL) (3 ' 2) 
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Equation (3-1) can be integrated with respect to y in the interval 0 < y ^ h to 
yield 


or 


U h . 3 aT ,3 . 9 yU 

. , o h 3p x ,3T\ f h ,dp N ^ , o 

pC ' 127 ( * } = [ T27 ( 3^ + “h~ 


(3-3) 


(il) --L 

v 3x' pc 


12y 3x 


3 2 ° 2 

,3p.2 y_o 

''Iv' ' U 


U h 
o 


h 3 3p 


12y 3x 


(3-4) 


Normalizing and utilizing the relationship 


(1/y) (3p/3?) = (h - K)/(h 3 ) 


the following can be obtained 

3T/3? = [y(C)/K] { [ (h - K) 2 /(h 3 )] + [ 1/ (3h) ] } (3-5) 

where 

6y U L(T - T ) 

_ O O O 

3-1-2 Lubricant Flow 

The flow pattern at the interface of a pumping ring was somewhat problematic. 
The nature of this flow can' be visualized in Figure 3-1. As shown in Figure 
3-la, at zero or low upstream pressure, the flow is mainly forward, with only a 
small pocket of fluid circulating near the inlet of the film. This recirculat- 
ing flow is induced by the adverse hydrodynamic pressure gradient prevailing 
near x = 0. When the level of p^ rises, more and more of the forward flow near 
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the stationary surface is blocked and eventually reversed (Figure 3-lc) by the 
static pressure gradient induced by p^. Eventually, when the upstream pressure 
is sufficiently high, fluid from the sealed chamber begins to leak backwards 
along the stationary surface. When the magnitude of p^ becomes very high (Fig- 
ure 3 -Id), most of the flow is backward, with only thin, vanishing layers of 
forward and recirculating flows maintained near the moving surface. 

Thus, in general, there are three layers of flow possessing the following char- 
acteristics: 

• Forward Flow - This flow is along the moving surface. It enters at a 

temperature, T , that prevails at x = 0 and is heated on its way to the 

reservoir to some T at x = L. 

max 

• Recirculating Flow - This flow also enters at a temperature, T , but 
various portions of that flow penetrate only part way into the film 
before they are reversed and returned to their source. It should be 
noted that the bulk of the flow recirculates near the entrance, where h 
is large. It undergoes relatively little viscous shearing, resulting in 
low energy dissipation to the fluid. 

• Reverse Flow - This flow originates in the reservoir entering at a 
temperature, T^, and is heated while traveling upstream. Its maximum 
temperature is reached at the entrance to the pumping ring, at x = 0. 

In view of these characteristics, only fluid which transverses the whole length, 
L, i.e., only the forward and reverse flows, was considered instrumental in 
carrying away the dissipated heat. Since the bulk of the intermediate layer 
recirculates near the entrance where the temperature differential is relatively 
small, its effect is left out of the heat balance. This treatment represents a 
conservative approach because inclusion of the recirculating flow would yield 
lower temperatures and thus safer operating conditions than those predicted by 
the present method. 
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The first task was to find expressions for the three flow regimes in terms of 
ring geometry and its operating conditions. Defining a dimensionless transverse 
coordinate and a dimensionless velocity by 

u(y) = (u/U q ); n(?) = [y/h] 

the velocity from Equation (3-2) may be expressed in dimensionless form as 

u =. 3(1 - K/h) T\Cn - 1) + (1 - n) (3-6) 

The line of zero velocity, or the line below which all fluid flows forward and 
above which the flow is backward, is easily obtained from Equation (3-6) by 
writing u = 0 which yields the locus 

n(Ol - n*(?) O 7 ) 

1 u=*0 3[ l-K/h( £) ] 

At a given ?, the total flow contained between the moving surface and any point 
(B, ?) line is given by 

’/'(t 1 >?) = h ^ u(n ’ ) dn ’ 

and after integration yields 

*0i,$) = hn (H - l) 2 + Kti 2 (3/2 - n) (3-8) 

By assigning to <Kb, 5) various constant values, the flow streamlines are 
obtained. The net flow, of course, is contained between n = 0 and t| = 1 or 

*d,5) = = (K/2) 

which is shear flow at h = K where (dp/d?) = 0. 

The film thickness is given by 
h = h 2 [l + (a - 1)(1 - ?)] 
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So that the value of (K/2) for an isoviscous, linear slider is given by 


K 2a 2, 

r"rri- a - ap f h 2 > 

h 2 


where 


a = 


h 2 ' h 2 


, A 


12 / 


For the variable viscosity case where y = y(£)> the value of K is 

■f 


K = 


. . Udg 

P f + J ~ 2 

O*' h 


s $ 


(3-9) 


(3-10) 


The individual streamlines are obtained by assigning different constants to 

<Kti, 5) in Equation (3-7). A sample plot of such streamlines for a = 3 and 

2 

different values of p^h^ is shown in Figure 3-2. Reverse flow will commence 
when the dividing streamline between the forward and recirculating flow reaches 
£ = 1. This, from Equation (3-7), occurs at . 

* 1 

n (1) = — “ = 1 (3-11) 

3[ 1 - K/h 2 ] 

~ 2 

or at (K/h^) = 2/3. At values of K/t^ < 2/3, the tangent point of the recircu- 
lating envelope at £ = 1 will lie below tj* - 1, opening up a passage for reverse 
flow. Below the recirculating envelope the fluid will, as shown in Figure 3-3, 
flow forward at a rate of 


q p = <Mti*(1); 1] 


whereas, above the envelope there will be reverse flow at a rate of 


q R = q p - q N£T = ^[Tl*(l), 1] - (K/2) 


(3-12) 
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3.1.3 Modeling of Thermal Problem 


As stated previously, the viscous heating is assumed to be convected away only 
by the fluid that passes the interspace, be it to the right or to the left. For 
this purpose, the interspace is modeled so as to be filled with these two 
streams, i.e. a forward flow, q_, and a reverse flow, q n , as shown in Figure 3-4. 
By designating 

f=q F /(q p +q R ) (3*13) 

f0(x) will be the heat convected to the right by q^, and (1-f) 6 the heat 

r 

connected to the left by q D . Likewise, temperatures will be averaged across the 

film (i.e., in the y but not in the x direction) with T_(x) and T (x) represent- 

r K 

ing the temperature profiles generated in the two flow layers, q F and q R . The 
overall average temperature profile will be obtained from 

T(x) = f T p (x) + (1 - f) T r (x) (3-14) 

The viscosity at each x station will be based on this averaged temperature, i.e. 

V(x) = li[T(x)]. 

The expression for the viscous heating (Section 3.1.1) is exact in the sense 
that the losses are calculated over the exact velocity profile, including the 
region of recirculating flow, namely- 

<Kx) = y(x) Q / h (du/dy) 2 dy (3-15) 

where u (x,y) is the velocity profile shown in Figure 3-3. Consistent with the 
flow model, y(x) in the calculation of this viscous shear will be that corre- 
sponding to the average temperature T(x). 

It should be noted that, whereas with no reverse flow (q_ = 0), the temperature 
at the inlet to the film is a constant equal to T , this is no longer true when 
reverse flow sets in. The reverse flow, starting at an initial temperature, T^, 
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will reach a maximum at x = 0 and, since temperatures are averaged across y, the 
inlet temperature will be some T ' > T . 

3.1.4 Calculation Procedure 


The procedure to be followed in calculating thermal effects in the fluid film 
is, in view of the previous discussion, as follows: 


a. The differential equation to be solved is 




h - K 


where due to the variable viscosity y(5) 


h 



b. If (K/h^) > 2/3, there is no reverse flow and 


dT _ u(Q 
d£ ~ K 



with T = T at 5 = 0. 
o 


(3-16) 
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c. If (K/h^) < 2/3, there is reverse flow, and the following quantities 
have to be established in order to account for both the forward and 
backward streams : 


k 

n (1) 



K p = 2h 2 n*(n* - l) 2 + 2K(n*) 2 [3/2 - n*] = 2q p 


Kr = Kp - K = 2q R 

Kl = Kp + k r = 2Kp - K = 2(q F + q R ) = 2q T 

dT = p(Q 
dC ^ 


(h - K) 2 _1_ 

~3 ' 

h 3h 


(3-17) 


Equation (3-17) needs a solution with two different boundary conditions: 

• For Tp (5) , T = T q at \ = 0 

• For T r ($), T = T f at $ * 1 


The averaged temperature, any given is thus 


T(€) + 


t r k f + 




(3-18) 


where 

pc C 2 (T - T ) 

T = 6u U L 
o o 

A sample solution for the temperature profile in the fluid film for various 
values of upstream pressure is given in Figure 3-5. The solutions are for the 
case of equal upstream and downstream boundary temperatures . Curves which start 
at T = 86°F are those without reverse flow and thus have what may be called a 
conventional profile. However, at p^ > 1.5, there is reverse flow and, due to 
the averaging of temperatures across y, the inlet temperature is higher than 
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86°F. It is interesting to note that, due to the cooling effect of the backward 

flow, the maximum temperatures in the film actually decrease as the back flow 

increases above a certain level. Thus,.,. for p.,. = 0, . T„ . = ,10,7°F,. and for 

p^ = 3.0, = 92°F and occurs not at the trailing but 'at the leading edge of 

the film. The highest temperatures occur at some intermediate combination of 

forward and backward flows; in this particular example.. it happens at p^ = 2.0 

with T reaching 120°F. 
max & 


3.2 Variable Velocity and Squeeze-Film Effects 


Throughout the previous discussions,, the rod velocity was considered to be 
constant, given by U q = 2sf. This, of course, represents the average velocity 
over each half cycle. In actuality, the rod, driving a crankshaft, moves with a 
variable velocity given by 


u = u cosir2ft 
max 


where 


u 

max 


fs = ?rU /2 
o 




Consequently, a normal velocity component and squeeze film forces are imposed on 
the ring, as shown in Figure 3-6. These are generated by a variation of the 
hydrodynamic forces and of film thickness, as a function of the variable veloci- 
ty. The normal velocity introduces perturbations on nearly all the relevant 
quantities, such as film thickness, pressures, flows, extent of cavitation, and 
others. « •• 


When variable velocity and squeeze-film- effects ; j are included, the Reynolds 
Equation becomes: 


(3/30 [h 3 (8p/80] = (tt/ 2) [cosf(3h/90 +-a(3h/-3t) ] 


(3-19) 
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where o is the squeeze-film parameter given by a = 4L/s, and T is the dimension- 
less time t = 2nft. 


The relevant boundary conditions are: 


• For no caviation: p(0) = 0; p(l) = p^ 


= 0 

The one-dimensional, transient Reynolds Equation, given by Equation (3-19), may 
be solved with the use of the implicit time-transient method. 3h/dt may be 
written as: 


For cavitation: p(l) = p^, 


> 

(3h/ 3 t ) ; ^ h (i) ” h (i _ 1 )^ At 


where the subscript i in parentheses denotes the i t ^ 1 time step, and At is the 
interval between the i t ^ 1 and (i-l)*"* 1 time steps. Hence 


(i) h 2(i) 


+ Ah 


(i) 


(1 - 5 ) 


and the discretized form of Equation (3-19) is written as 




(i) 


- h 
At 


(i - 1) 


(3-20) 


If it is assumed that h^_^ is known, then Equation (3 _ 20) can be integrated 

’ w th 

analytically to obtain as a function of £ at the i time step. The elas- 

ticity equation, Equation (2-6), remains unchanged except for the time depend- 
ence of p, and may be coupled with the solution to Equation (3-20) and solved at 
each time step as described previously for the steady-state solution. Values of 
h 2 and Ah at the previous time step (^(i-l)’ ^(i-1)^ were usec * i n evaluating 
h^_^ in Equation (3-20) and for initial guesses for h^^ and Ah^^ for use in 
the secant method. 


40 



Quasi-static solutions (o = 0) were obtained at the middle of the forward stroke 
(t .= 0) to start the marching procedure. Solutions were then computed at 
successive time steps until periodic solutions were obtained. 

Sample solutions were run on a pumping ring having the following dimensionless 
characteristics : 

a = 0.0282, S = 0.2569, e = 0.518, L = 1.567, a = 4.0 

Three cases were considered, at different combinations of and p Q ; the 

P q = 4.0 solution also represents a case where the ring is clamped on the back- 
stroke. The three cases were all solved for steady-state (u = u q ), quasi-static 
(o = 0), and transient conditions to show the effects of squeeze film on ring 
performance. The value of a - 4.0 is quite large for practical applications as 
it corresponds to a stroke, s, equal to the bearing land, L. For present appli- 
cations, a is generally less than 1. The large value was used to exaggerate the 
effects of squeeze film for purposes of illustration and interpretation. Graphs 
for film thickness, flow, and pressure profiles are given in Figures 3-7 through 
3-14, whereas the values of the component flows are itemized in Table 3-1. Note 
should be taken that t = 0 corresponds in these plots to rod position at u = 

u , i.e., at the midpoint of its forward stroke, 
max r 

Considering first the effects of variable velocity alone (quasi-static 
solution), the plots show the following trends: 

• The film thickness variation during the forward stroke follows the sinu- 
soidal shape of the velocity curve. During the backstroke velocity, 
variation has no effect whenever cavitation commences at the trailing 
edge and a small effect when cavitation is located at £ < 1 (Figure 3-9). 

• The flow curve follows the shape of the velocity curve throughout the 
cycle, except, of course, when, due to a clamped ring, the flow during 
the backstroke is zero. In the particular case of Figure 3-14, the flow, 
after peaking at u max > became negative at the end of the forward stroke 
(u •* 0). However, upon commencing the backstroke, the ring clamped and 
flow ceased. 
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Fig. 3-12 Pressure Profiles as Function of Time for Transient Analysis 




Minimum Film Thickness 



Fig. 3-13 Film Thickness B 2 





I I I 1— — i t 1 1 I 1 

240 260 280 300 320 340 360 


Time 






TABLE 3-1 


EFFECTS OF VARIABLE VELOCITY AND SQUEEZE FILM ON FLOW 

a = 0.0282 8 = 0.2569 e = 0.518 L x = 1.5614 

a = 4.0 


p 

" o 

p f 

JL 

Analysis 

FORWARD 

REVERSE 

NET 

Q 

AQ, % 

Q 

AO, % 

Q 

AQ, % 

3.0 

0 

S 

0.4542 

Ref. 

0.1816 

Ref. 

0.1363 

Ref 



QS 

0.4687 

+ 3.2 

0.1811 

0.28 

0.1438 

+ 5.5 



T 

0.4231 

- 6.9 

0.1731 

4.7 

0.1250 

- 8.3 

3.0 

1.0 

S 

0. 3392 

Ref . 

0.4080 

Ref. 

- 0.0344 

Ref. 



QS 

0. 3524 

+ 3.9 

0.4059 

- 0.8 

- 0.0268 

1 22 | 



T 

0.3020 

- 11.0 

0.4110 

+ 0.7 

- 0.0545 

• 1 58 | 

4.0 

2.4 

S 

0.0948 

Ref . 

0 


0.0474 

Ref. 



OS 

0.1348 

+ 42 

0 

— 

0.0677 

+ 43 



T 

0.0711 

+ 25 

0.1165 

— 

0.0227 

- 148 


* S - Static 


QS - Quasi-static 


T - Transient 





































• Variable velocity has the effect of increasing the absolute values of all 
flows; that is, even though it boosts both the forward and back flows, 
the net flow is increased. Table 3-1 shows this increase to be of the 
order of several percent for the case of p^ = 0. The high percentage 
changes shown for other cases should not be given too much signficance 
since they center about very low flow levels when a small variation is 
apt to produce deceptively large effects in terms of percentages. 

• The pressures are, of course, positive over the entire ring during the 
forward stroke, as shown by the 0 - tt range of the constant lines on Figure 
3-11. Their magnitude, too, follows the velocity curve. During the 
backstroke, tt/2 < t < 3ir/2, cavitation sets in as high negative veloci- 
ties are approached (t + rr); they tend, however, to disappear at the 
beginning and end of the backstroke, when the negative velocities are 
low. 

The inclusion of squeeze-film effects have the following effect: 

• The largest film thickness does not occur at u . The peak is delayed 

so that it is reached after u . Likewise, the onset of constant h 

max 

during the reverse stroke does not commence at u = 0, but is also 
delayed. Thus, squeeze films have the effect of producing a phase shift 
in the film thickness curve. The shift is fairly large, varying from 50° 
to 90° in the three cases considered. However, the absolute values of h 
seem not to be affected. 

• The effect on the flow curves is similar in that there is a phase shift 
with little change in their magnitude. Both the forward and net flows 
are reduced in comparison to either the steady-state or quasi-static 
solutions; the amount of back flow tends to increase under the influence 
of squeeze film forces. Since the quasi-static solution gives higher 
flows and the static solution gives lower flow than the constant parame- 
ter approach, the inclusion of variable velocity without squeeze film 
effects would yield errors of at least 10%. 
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•Unlike the behavior in the quasi-static solution, the pressure curves 
are not identical over the O-TT and tt- 2 halves of the cycle. The pres- 
sures are lower when squeeze film effects are included, except over the 
latter half of the backstroke, when a closing gap helps to increase the 
pressures above those of the quasi-static analysis. 

3.3 Starvation 

An estimate of the effect of starvation can be obtained by first considering 
the clamped case where the area under the bearing land, 0 < x < L, is 
predicted to be completely dry at the end of the backstroke. The average film 
thickness in this area is Ah/2, and the average Couette flow per unit of 
circumferential length is U Q Ah/4. The time, t s , for the starved volume to 
fill up would thus be the cross-sectional area AhL/2 divided by the flow per 
unit of circumferential length. Thus, 

t g = (AhL/2) /(U Q Ah/4) 

The distance traveled prior to flooding is U 0 t s = 2L. This would reduce the 
effective forward stroke by an amount equal to 2L. The value of t s given above 
will be somewhat low in that it does not account for the development of a 
resisting pressure gradient that will occur as the starved volume fills up nor 
does it account for any inertial effects in the entrance region that could 
impede the start of the filling process. In order to account for these 
effects when analyzing experimental data, multiply t s by a factor X which will 
in general be greater than or equal to 1. Thus, 

s e ff = s - 2LX 

and the dimensionless flow rate K e ff becomes 

K eff = K ( s e f f ^ = K [i " X(2L/s)] (3-21) 

The starvation process is, of course, much more complex than that treated 
here, especially in the unclamped case where there is a partial film and the 
film is exposed to the sealed pressure, p Q , so that the cavity could be 
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trapped in the middle. In order to understand more fully the influence of 
starvation on pumping ring performance, transparent model experiments should 
be performed so that the size and behavior of the partial film can be observed 
and modeled. In the interim, a simple model was adopted to estimate the 
effect of starvation on pumping for either the clamped or cavitating, 
none lamped case. 

If the film cavitates during the backstroke at £ = £ c > the film thickness at 
the point, h c , is 

h c = h 2 + Ah(l - £ c ) 

and the void volume (volume of gas or vapor in the cavitated region) per unit 
of circumferential length, V c , would be 

V c = (h 2 + Ah - h c ) L £ c /2 = L C c * Ah/2 

If it is assumed that the void volume must fill up before significant pumping 
can begin and if it fills based on Couette flow at the average film thickness 
of the cavitated region, the incoming flow per unit of circumferential length 
would be U Q (h 2 + Ah + h c )/4, and the starvation time, t s , would be V c divided 
by that flow 

t s = (2U c 2 Ah)/[U Q (h 2 + Ah + h c )] 


As for the clamped case, the effective forward stroke would be reduced by Xu Q 
t s . If the dimensionless flow in the forward stroke is denoted by Kf and the 
backstroke by K^, then 


K = K (s - XU t )/s 

t err f os 


= K r 


= K. 


1 - 2 \ 


1 - 2X 


Hr) 

0) 


r 2.. "l 

£ Ah 


h_ + Ah + h 
2 c 

g 2 Ah 

2h 0 + Ah (2 - ^ ) 
^ c 
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and 


^eff = Kf e ff Kg (3-22) 

For the clamped condition, Kg =0, 5 C = 1 and h2 = 0; thus Equation (3-22) 
reduces to Equation (3-21). 

In order to show the effects of starvation, the data for a clamped carbon 
graphite ring is given in Figure 3-15. The dotted curves include the effects 
of starvation for a value of X = 1. It can be seen that, although the inclu- 
sion of the effects of starvation reduces the discrepancy between theory and 
experiment, it does not resolve the lack of agreement. The effect of values 
of X > 1 will be shown later when comparing theory to experiment. 

As a result of the uncertainty regarding the mechanisms of the starvation 
process, the effects of starvation are not included in the parametric studies 
in Section 4.0. It is recommended, however, that Equation (3-22) be used to 
estimate potential effects of starvation when designing pumping rings. The 
dimensionless starved net flow is given as an output to the computer program 
RING listed in Appendix B of this report. 

3.4 Nonparallel Contours 


Essentially, the contours considered are tapered surfaces with a constant 
slope as shown in Figure 3-16. The reference clearance is that at the trail- 
ing edge of the ring; the leading edge clearance is designated by Cm- The 
slope parameter, 6, is then given by 

6 = (dh/dx) = (C - C m )/L (3-23) 

In terms of the nondimens ional quantities h and £, the slope is given by 

6 = (dh/dC) = (C - O/C = (1 - CJ (3-24) 

M M * 


54 



12 


Ol 


11 h 
10 


1 

0 



Fig. 3-15 Effect of Starvation on Performance of Carbon Pumping Ring 
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Fig. 3-16 Geometric Taper on Pumping Ring 
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The clearance ratio Cm = (Cm/C) in terms of this 6 is then given by: 

C M = 1 - L6/C (3-25) 

The effect of using tapers on the performance of pumping rings is discussed in 

Section 5.0. Briefly, it can be said that the use of tapers has little effect, 

except in cases of low clamping load, p , or high values of E, i.e., when there 

o 

is little elastic deflection of the ring. Clearly, when there is no 
deflection at all, a taper would constitute the sole mechanism of generating 
hydrodynamic pressures. In practical applications where deflections are 
present, the effect of tapers is minimal. 
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4.0 PARAMETRIC STUDY 


A study was made regarding the impact of various structural and operational 
parameters on the performance of the pumping rings. These parameters include: 

• Modulus of Elasticity (E) 

• Poisson's Ratio (v) 

• Clearance (C) 

• Length (L & Lj) 

• Ring Thickness (t) 

• Taper ( 6 ) 

• Loading (p Q and e). 

The computer runs that provide the results of the parametric study are based 
on the constant-parameter approach without starvation but including the 
effects of the backstroke and attendant cavitation. Thermal and transient 
effects were omitted because, as shown on Figure 4-1, their effect is not of 
the sort as to qualitatively change the behavior of the pumping ring. The 
constant parameter analysis should be sufficient to reveal the essential 
features of the parametric relationships without the undue complications of 
the more elaborate analysis. 

Since it would ultimately be desirable to optimize the design of pumping 
rings, the question arises as to what constitutes an optimum. Given the 
purposes for which these devices are customarily used, an optimum ring was 
deemed to be one which, within given constraints such as reasonable values of 
p o , C, etc., can maintain the highest possible sealing pressure without exces- 
sive wear. This implies the highest values of p^ and a low, preferably zero, 
frictional force F. Thus, a ring design which generates high p^ and just 
clamps shut upon reversing the stroke would fulfill this criterion. An addi- 
tional quantity of interest, perhaps, is also the maximum flow Qq (flow at p£ 
= 0), that such a ring can produce. Thus, in the plots that follow, the items 
of pf , F, and Q q will be the items against which ring performance will be meas- 
ured. 

The parametric study was conducted by first establishing a standard or refer- 
ence design and then by varying its parameters, one at a time, to values below 
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and above the reference quantities. The parameters for this standard ring are 
given in Table 4-1 along with selected departures from the standard set. 
Including the standard, most parameters have four computed points from which 
optimization plots can be made. The results of the parametric study are 
summarized in Table 4-2; more relevant plots are portrayed in Figures 4-2 
through 4-11. In Figure 4-4, the reason the film thickness h 2 drops with an 
increase in clearance, is due to the drop in Pf m » a high p^ usually helps to 
contract the clamping pressure and thus maintain a high h2* In Figure 4-6, 
the explanation for an increase in flow with a rise in length L is that with a 
greater L the hydrodynamic effects are increased, and thus, also the flow. As 
seen, four parameters; the lengths L and Lj, taper 6, and Poisson's ratio have 
little effect on performance. Of the two important variables, namely, clear- 
ance and clamping force, the first has to be raised above its optimum for 
maximum p£ , and the latter reduced below the optimum, if large frictional 
forces are to be avoided. Since the clamping force is given by the product of 
p 0 and e, Figure 4-12 shows the effects of using the same clamping force, 
31.96 KN/m (182.5 lb per in.) of circumference, but by varying the relative 
values of p Q and e. As seen, higher values of p^ are achieved when p Q is high, 
while higher values of are achieved when e is high. The differences, 
however, are not striking. 

It should be pointed out that the particular optima recorded in Table 4-2, 
such as t = 2.54 mm (0.1 in.), C = 25.4 microns (1 mil), etc., are valid in the 
range of parameters characterizing the particular ring specified in Table 4-1. 
Thus, for example, for a material with a much higher E value, optimum results 
would be achieved with lower values of clearance and ring thickness. The 
opposite would be true for a material with a value of E much below the 34.5 GPa 
(5 x 10^ psi) assigned to the standard design. 



TABLE 4-1 


STANDARD DESIGN FOR PARAMETRIC STUDY 



e 


T 

E 




s 

f 

R 

L 

L. 

C 

t 

6 


10.3 MPa (1500 psi) 

3.175 mm (0.125 in.) 

49°C (120°F), ii = 59 x 10 ^ Pa-sec, (8.6 x 10 ^ Reyns) 
3.45-10 4 MPa (5-10 6 psi) 

0.36 

12.7 mm (0.5 in.) 

35 Hz 

12.7 mm (0.5 in.) 

6.35 mm (0.25 in. ) 

10.2 mm (0.4 in.) 

0.019 mm (0.75 x 10 ^ in.) 

1. 9 mm (0. 075 in.) 

0 

VARIATIONS IN PARAMETERS 


-3 -3 -3 

C: 0.0063 mm, 0.0127 mm, 0.0381 ram (0.25 x 10 in., 0.5 x 10 in., 1.5 x 10 in.) 

E: 6.895- 10 3 MPa, 20.7-10 3 MPa, 68 . 95 -10 3 MPa , 207 -10 3 MPa (10 6 psi, 3 -10 6 psi, 

10*10^ psi, 30*10^ psi) 

L: 5.1 mm, 7.6 mm, 10.2 mm (0.2 in., 0.3 in., 0.4 in.) 

P q : 5.17 MPa, 6.895 MPa, 13.79 MPa (750 psi, 1,000 psi, 2,000 psi) 
v: 0.25, 0.5 

t: 1.27 mm, 2.52 mm, 3.81 mm (0.05 in. * 0.1 in., 0.15 in.) 

L 6.35 mm, 15.2 mm (0.25 in, 0.6 in.) 

6: -10" 3 , -2- 10~ 3 , -3- 10~ 3 
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TABLE 4-2 


OPTIMUM PARAMETERS FOR STANDARD RING 
(See Table 4-1) 


ITEM 

RANGE 

EVALUATED 

For Pfm 

FOR Q 

o 

SUBJECT TO F=0 

E 

psi 10 ^ 

1 - 30 

5.5 

6 


MPa 10~ 3 

6.89 - 207 

41.3 

41.3 

n 



No Effect 

No Effect 

c 


' 6:25 - 1.5 

1.0 

1.0 



6.35 - 38. £ 

25.4 

25.4 

a 

In 

0.2 - 0.4 

No Effect 

Highest 

B 

mm 

5.08 - 10.16 



B 

in 

0.25 - 0.6 

No Effect 

No Effect 

B 

mm 

6.35 - 15.2 



B 

in 

I 

0.100 

0.102 

B 

mm 


2.54 

2.64 

6 


0 - (-3-10' 3 ) 

No Effect 

0 

P o* 

psi 

750 - 2,000 

1100 

1050 


MPa , 

5.17 - 13.8 

7.5 

7.20 


A ’ 





, ,*See Also„ Fig . 4-12. 
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Fig. 4-2 Effect o f E on Performance Standard King (See Table 3- 1 ) 
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5.0 EXPERIMENTAL PROGRAM 


Following the optimization of the candidate rings, tests were run on three 
rings to verify their performance and their agreement with theoretical pre- 
dictions. The test rig and the experimental procedures were the same as those 
used on the tests described in Ref. [l]. 

The pumping ring program, given in Appendix B based on the analysis performed 
here, has been used to select ring designs for experimental testing. Initial- 
ly, a set of studies was performed with steel, babbitt, and Rulon J rings as 
examples of a high, medium, and low modulus pumping ring. Overall dimensions 
were selected consistent with the experimental test rig. The principal design 

variables to be evaluated within the applied pressure, p , the clearance, c, 

o 

and the average thickness, t. A target value of p of 1500 psi was used; 

o 

however, it was found not feasible to design to this pressure with the low 
modulus Rulon ring. Sample geometries that have been arrived at are shown in 
Table 5-1. The general configuration of the pumping rings submitted for test 
is shown in Figure 5-1. 

The principal criterion used involves being able to significantly deflect the 
ring to obtain pumping action without excessive clamping during the back- 
stroke, which would result in excessive wear. 

Based on the criterion above, the design of the steel ring would require rela- 
tively small thickness and a small clearance to obtain suitable compliance and 
deflection characteristics. The tolerances required seem excessive for the 
present application, and it is thought that steel rings would only be applica- 
ble in situations requiring much higher pressures. It was thus decided to 
confine testing to medium and low modulus materials ranging between babbitt, 
carbon graphite, and Rulon. The test rig and experimental procedures were the 
same as those used on the test. 

One problem area revealed in the course of testing concerns the clamping load. 
Difficulties were encountered in testing the carbon and Rulon rings, and often 
such attempts at testing led to breakage of the specimens. It was also noted 
that some of the babbitt rings were severely worn, not at the downstream end 
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TABLE 5-1 


SAMPLE GEOMETRIES AND CONDITIONS FOR 
PARAMETRIC STUDY OF CANDIDATE TEST RINGS 


L = 6.73 mm ( 0.265 in.), L^ = 7.62 mm (0.30 in.), 6=0, s = 38.1 mm (1.5 in.), 

f = 35 Hz. 


ITEM 

STEEL 

BABBITT 

RULON J 

R, mm 

9.5 

9.5 

9.5 

(in. ) 

(0.38) 

(0.38) 

(0.38) 

E, MPa 

co 

o 

H 

r- 

O 

cs 

51- 10 3 

1.72* 10 3 

(psi) 

(30 x 10 6 ) 

(7.5 x 10 6 ) 

(0.25 x 10 6 ) 

V 

0.3 

0.36 

- 0.46 

p , MPa 
o 

10.34 

10.34 

5.17 

(psi) 

(1500) 

(1500) 

(750) 

C, microns 

8.96 

19 

50.8 

(mils) 

(0.35) 

(0.75) 

(2) 

t, mm 

1.27 

1.90 

3.81 

(in. ) 

(0.05) 

(0.075) 

(0.15) 
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where wear might be anticipated but far upstream, ahead of the actual pumping 

area, as shown in Figure 5-2a. This Led to the conclusion shown in Figure 5-2b 

that the high-pressure oil used for clamping leaked past the O-ring and, by 

pressing against bracket A, effectively produced a sealed chamber over the OD 

of the ring. The pumping ring was thus loaded with a pressure, p , over the 

o 

entire length (L^ + L) instead of only the narrow region, e, covered by the 
O-ring. This unpresidented pressure loading gave the ring the flared shape 
shown in Figure 5-2b which led to high wear upstream of the active pumping 
area and to breakage of the weaker rings. 

The problem was resolved by opening relief passages behind bracket A so that 
any fluid leaking past the O-ring would be scavenged outside. Since the leak- 
age past the O-ring occurs only when p. approaches p , such leakage should not 

r o 

affect the build-up of sealed pressure p^ at levels below p^. With the relief 
grooves in place, all unsatisfactory tests were repeated; no further difficul- 
ties in testing the Rulon and carbon rings were encountered. It should also 
be noted here that the sealed pressure was not observed to significantly 
exceed the loading pressure p Q . This is believed to be a result of leakage 
past the O-ring, which occurs when p^ starts to exceed p^, as indicated in 
Figure 5-3. 

The rings tested were made of either babbitt, carbon graphite, or Rulon. The 
geometry and dimensions of the various rings tested are as summarized in Table 
5-2 and, except for one item (Item VII), had the following two dimensions in 
common: 

Shaft Diameter - 19.05 mm 

Inside Diameter of Back Section (over length L^) - 19.3 mm 
The viscosities of the oil used are those given in Figure C-l of Appendix C. 

5.1 Tests with Large Babbitt Ring 

The 19.05-mm diameter babbitt ring was tested at three frequencies, 60, 35, 
and 10 Hz, and at two strokes, 50.8 mm and 25.4 mm. The parameters investi- 
gated were the pumping ring length (L), clearance (C), and the effect of a 
taper (6). Tables C-l through C-4 in Appendix C give the detailed results of 
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High Pressure Oil 




Fig. 5-2 Leakage of High Pressure Oil Past 0-Ring 
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Fig.. 5-3 Schematic of Leakage Past 0-Ring at 
High Sealed Pressures 


TABLE 5-2 


LIST OF PUMPING RINGS TESTED 


oo 

o 






(.843) 


(.267) 

V 

Carbon 

9.52 

(.375) 

22.07 

(.869) 

7.56 

(0.298) 

6.78 

(.267) 

VI 

Rulon 

9.52 

(.375) 

24.15 

(.951) 

7.56 

(0.298) 

6.78 

(.267) 

VII 

Babbitt 

t 

6.00 

(.236) 

13.92 

(.548) 

5.59 

(0.22) 

5.08 

(.200) 


( 1 . 0 )* 


( 0 . 8 ) 

42.55 

(1,675)* 


8.89 

(0.35) 


3.75 


* Average of two rings. 











Che tests. The self-pressurization runs represent an arrangement in which the 
clamping pressure was provided by the sealed pressures generated by the pump- 
ing ring. This represents the case of p = p at any instance. It turned out 

o f 

that, in order to start the untapered pumping ring, a priming pressure was 
required (that is, when at the start, p^ = 0, no self-clamping was possible). 
Some initial p^ Q > 0 had to be provided by external means in order to enable 
the pumping ring to apply self-pressurization and build up a p^ beyond the 
initially supplied priming pressure, p^ Q . The tapered rings were found to be 
self-priming. 

With regard to the effect of the several variables tested, the plots suggest 
the following: 

• Effect of Clearance - Table 5-3 shows the changes in maximum flow and 
pressure induced by changing the clearance from 11.4 microns to 19 
microns. These changes are minimal, indicating that the optimum lies 
within the range of values. 

• Effect of Length - As shown in Table 5-4, there is a definite advantage 

to a higher length, L, with regard to maximum flow rates; but the 

shorter rings produced higher maximum levels of sealed pressure, p . 

rm 

• Effect of Taper - Table 5-5 confirms what can be intuitively deduced as 
the desirability of having a geometric taper. At high clamping pres- 
sures when ring deflections are large, the geometric taper is counter- 
productive. However, at p^ = 3.45 MPa, a taper produces higher flow 
rates and higher levels of maximum pressure. 

5.2 Tests with the Rulon and Carbon Ring 

The detailed results of the tests on the Rulon and carbon rings are given in 
Tables C-5 and C-6 in Appendix C. In both rings, the higher clamping pres- 
sures produced higher reservoir pressures, p , but the lower values of p 

fm o 

yielded higher maximum flow rates. This trend was already noticeable with the 
babbitt ring, but as the value of E drops, it becomes more pronounced so that, 
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Q 0 , gra/min 


MPa 


Q 0 , gm/min 

Pfm’ m * 

= 0. 0114 mm 

C = 0.019 mm 

C = 0.0114 mm 

C = 0.019 nun 

(0.45 mil) 

(0.75 mil) 

(0.45 mil) 

(0.75 mil) 

43 

42 

9 

9 

21 

22 

9 

9 

4 

4 

8.8 

8.8 

50 

51 

7.2 

7.2 

24 

27 

7 .'.2 

7.2 

4 

4 

7.2 

7.2 

7 

7.5 

8.8 

9.2 
























TABLE 5-4 


EFFECT OF REDUCED LENGTH L 
RUNS I AND III 


P o 

MPa (psi) 

s 

mm (In. ) 

f. 

Hz 


Q 0 , gm/min 

p f„' 

L= 6. 78 mm 
(0.267 in.) 

L=4 . 83 mm 
(0.19 in.) 

L=6. 78 mm 
(0.267 in.) 

L=4 .83 mm. 
(0.19 in.) 




■ 





6.895 

50.8 

60 


50 

34 

mm 

7.8 

(1000) 

(2) 





mm 




35 


25 

17 

E9 

7.8 



...IQ 


.4 

2.5 

WSM 

7.0 

3.45 

50.8 

60 


59 

51 

3.8 

mm 

(500) 

(2) 






■ ■ 



35 


28 

27 

3.8 




10 


3.5 

3.5 

3.8 

n 

8.62 

Bjp5iiS3| 

60 


70 

60 

8.8 

8.0 

(1250) 

mm 
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5-5 


T OF TAPER 
I AND IV 


V 

gjn/min 

P fm’ 

MPa 

= 0 

6=-3.75*10' 3 

6 = 0 

6=-3.75*10~ 3 





42 

32 



22 

17 

9.2 

9.2 

4 

4 

8. 9 

9.0 

50 

42 

9.2 

7.2 

24 

18 



4 

2 


7.5 

60 

o 

oo 

A 

3.8 

4.2 

28 

38 

3.8 

4.2 

4 

4 

3.8 

3.8 

7 

2.5 

8.8 

8.8 
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for the Rulon rings, Q at p = 1.72 MPa (250 psi) is two or three times the 

o o 

value of Q at p = 5.17 (750 psi). 
o r o 

5.3 Tests with the Small Babbitt Ring 


The results for the 6-mm radius babbitt ring are given in Table C-7 . The diam- 
eter was reduced roughly 1/3 as compared to the large rings. Neither the 

flows, Q , nor the maximum sealed pressure, p , were much affected by the 
o tm 

change in size. 
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6.0 DISCUSSION OF RESULTS 


Three sets of comparative results are presented in Appendix D. 

• The babbitt rings, which comprise the large size with two different 
clearances and a small size ring 

• The carbon graphite ring 

• The Rulon ring. 

The theoretical curves presented in Appendix D were obtained without the star- 
vation correction and without allowance for the fact that the sealed pressure 
does not exceed the loading pressure due to leakage past the O-ring. 

! 

A number of qualitative conclusions can be drawn from a study of the figures 
in Appendix D. The major discrepancy observable is the fact that the measured 
flows are substantially smaller than those predicted by theory at zero sealed 
pressure. This is believed to be partially due to effects of starvation as 
will be discussed later. For practical purposes, this discrepancy is not 
believed to be a major one in that the pumping rings will normally be used to 
develop a buffer pressure as opposed to being designed to deliver prescribed 
flow rate. The predicted maximum developable pressure at zero net flow was 
found to be in better agreement with the experimental data. Again, all of the 
curves that predict sealed pressure higher than the loading pressures should 
be truncated at the loading pressure value due to leakage past the secondary 
seal. 

It is important to note that all of the data presented in Appendix D was 
obtained with seals that were sized and designed with the analysis contained 
here. In all cases, as predicted, the seals did perform the task of develop- 
ing and maintaining the designed sealed pressures. Significant pressures 
could not be developed for steel rings where the theory indicated that suffi- 
cient deflection would not be obtainable to provide adequate pumping. 
Although the analysis presented so far does not accurately predict deliverable 
flow, it can be used reliably to size and design pumping rings. In order to 
obtain and improve fit to the flow data, a two-constant empirical relationship 
has been determined from the data for babbitt pumping rings to attempt to 
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account for starvation. The results of this fit when compared with exper- 
imental data are described below. 

6.1 Empirical Correction for Starvation 

An empirical relationship for the starvation factor, X, as a function of L/ s 
of the form 

X = 0.74 (s/L) 0 * 58 

has been obtained by fitting flow rates at zero sealed pressure for ring No. I 
at p Q = 8.62 MPa. Comparisons have been made for rings I - III over a range of 
frequencies, strokes, loading pressures, sealed pressures, and geometries. 
Values of X obtained from the above equations varied between 1.6 and 2.4 for 
all the cases shown in Figures 6-1 through 6-10. 

Comparisons of the data for a 50.8-mm stroke at 35 Hz in Figure 6-1 with those 
for a 25.4-mm stroke at 60 Hz shown in Figure 6-4 indicate that, even though 
the speeds (product of stroke and frequency) are fairly close together, the 
short stroke data in Figure 6-4 shows a factor of 3 less flow. This is very 
much in keeping with the starvation analysis that predicts that reduced flow 
will occur when the land length becomes an appreciable fraction of the stroke. 
If starvation is neglected, the theory would predict the flows to be solely a 
function of the product of frequency and stroke. Even though only a two- 
constant fit was used, the agreement between theory and experiment for all of 
the cases involving the small clearance babbitt ring looks reasonably good. 

Figure 6-6 shows the comparison between theory and experiment with the larger 
clearance babbitt ring. The shaded area shows the difference between pre- 
dictions for a 0.75-mil clearance and a 0.65-mil clearance. The difference 
being well within the limit of the measurements. In general, the data tend to 
indicate that the 0.65-mil clearance is probably closer to the truth. This 
appears to be particularly true when one looks at the data at 10 Hz. Again, 
the agreement is reasonably good and no additional constants were used in 
fitting the data for the large clearance ring. 
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Fig. 6-1 


Comparison Between 
(s = 50.8 mm; p Q = 


Corrected Theory and Experiment for Babbitt Pumping Ring I 
8.62 MPa) 
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Fig. 6 2 Comparison Between Corrected Theory and Experiment for Babbitt Pumping Ring I 
(s = 50.8 mm; p n = 6.89 MPa) 
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Fig. 6-3 Comparison Between Corrected Theory 
(s = 50.8 mm; p Q = 3.45 MPa) 
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Fig. 6-4 Comparison Between Corrected Theory and Experiment for Babbitt Pumping Ring I 
(s = 25.4 mm; p 0 = 8.62 MPa) 
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Fig. 6-5 


Comparison Between 
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Fig. 6-7 Comparison Between Corrected Theory and Experiment for Babbitt Pumping Ring II 
(s = 50. R mm; p q = 6.89 MPa) 











Fig. 6 8 Comparison tie tween Corrected Theory and Experiment for Babbitt Pumping Ring II 
(s = 25.4 mm; p Q = 8.62 MPa) 
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Figures 6-9 and 6-10 show comparisons between theory and experiment for the 
shorter length ring. The increased discrepancies indicate that the starvation 
correction depends more on the stroke than on the seal land. This could be an 
effect of inertia in the inlet region. In general, the starvation correction, 
along with the constraint that the end pressure doesn't exceed the loading 
pressure, provides much better agreement between theory and experiments but is 
still somewhat incomplete. 

6.2 Effects of Viscosity 

Since the viscosity cannot be measured at the interface (it is, in fact, meas- 
ured several inches upstream of the inlet) and since the clearances are very 
small, both parameters could, in actuality, differ from the assumed values. 
Major reductions in viscosity, such as that which would occur if the actual 
temperature were 100°F higher than the measured inlet temperature, would 
provide remarkably improved agreement between theory and experiment. This 
fact was observed in numerous parametric studies which are not presented here 
in that no justification has been found for this temperature rise. 

6.3 Suggested Design Procedure 

The most effective use of a pumping ring is believed to be for generating a 
buffer pressure to back up another seal. Thus, if one were to seal gas at 10 
MPa from oil at ambient pressure, the pumping ring should be used to develop 
the 10 MPa backup pressure in a buffer reservoir, thus alleviating the pres- 
sure gradient on the primary seal. This should provide prolonged life for the 
sealing system since the pumping ring is essentially a "light contact" device. 

Clearances should be selected based on shaft diameters and machining toler- 
ances. Values of C/R should be of the order of 10”^, thus resembling bearing 
clearances rather than the tight clearances (or interference) that one would 
generally find in a seal. 

The loading pressure should be of the order of the desired buffer pressure. 
Back leakage past the O-ring could thus provide a potential relief valve to 
protect against excessive pressure buildup in the buffer chamber. 
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The computer program RING (given in Appendix B) should then be used to arrive 
at a combination of material properties (E, v) and geometries (t, e, L, Lj) so 
that the computed value of Pf m is equal to the desired buffer pressure, and a 
small amount of clamping is predicted to occur during the back stroke. Allow- 
ances should be made for some wear to occur due to friction resulting from the 
radial shear force F. 

If the pumping ring is designed to generate its own loading pressure, p^ = p^, 
an initial taper should be used. Values of C^/C of the order of 2-3 have been 
shown to be adequate for obtaining self-pumping without the need for external 
priming. Caution should be used here in providing pressure relief to avoid 
overloading. 

Finally, if the pumping ring is used as a flow device, rather than a pressure 
generator, allowances should be made for the fact that the analysis tends to 
overpredict flow under many circumstances. In general, the "starved net flow" 
prediction generated by RING should be within a factor of . 2 of the flow that 
one would expect to obtain in practice, based upon the measurements reported 
here. 
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7.0 CONCLUSIONS 


A design analysis has been developed for predicting pressure-flow relation- 
ships for various input geometries, speeds, fluid viscosities, and elastic 
moduli of pumping rings. The analysis can be used to size and design pumping 
rings for various applications. Results of testing with Rulon, babbitt, and 
carbon-graphite rings at various loads, speeds, and geometries indicate that 
the analysis can consistently be used to design rings that work well over a 
fairly broad range of parameters. 

The design criteria for selecting the clearance would be dictated by tolerance 
and flow rate consideration. The applied loading, p Q , should be the order of 
the desired pumping pressure. The choice of materials and geometry are then 
dictated by elasticity considerations so that the deflected ring, under static 
load, clamps the shaft with a very small force to nearly eliminate the back 
.flow without resulting in excessive friction and ensuing wear and power loss. 

The greatest discrepancy between theory and experiment lies in the prediction 
of flow at zero sealed pressure, Q q . Predicted values of Q c are substantially 
higher than those observed experimentally in almost all cases. Predictions of 
the maximum sealed pressure at zero net flow, P^ m > differ from those observed 
experimentally, in that values of p^ substantially greater than the loading 
pressure, p Q , are frequently predicted but have never been observed. 

The analysis shows that starvation can have a significant influence on Q q , 
although not sufficient to completely explain the discrepancies, and that 
leakage past the O-ring could cause p not to exceed p Q . Further causes of 
discrepancy could be due to uncertainties in the clearance and the local 
temperature at the inlet to the pumping ring. 

It was found experimentally that the use of tapers on pumping rings enabled 
the rings to self-pump up to high sealed pressures. The untapered rings, in 
general, needed to be initially loaded until a sufficiently high sealed pres- 
sure could be generated to load the ring. No priming was necessary for the 
tapered rings. Theory indicates that performance of a tapered ring and an 
untapered ring at the same minimum film thickness are similar under signif- 
icant loading, but the tapered ring is predicted to pump even when unloaded. 
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APPENDIX A 


PRELIMINARY ANALYSIS OF PUMPING LENINGRADER SEAL 
'• •' NOMENCLATURE ' " 

6 Interference 



h 


x 


1 ’ 


X 


2 


P 


P 


g 


p i 


P 


o 


E 


v 


t 


R 

D = Et 3 /[12(l-v 2 )] 
f(x) 


y 

9 


Radial inward deflection 


Film thickness 
Coordinate variables 


Interface pressure 
Sealed gas pressure 
Interference pressure 
Ambient oil pressure 
Elastic modulus of seal 
Poisson's ratio 
Thickness of steel 


Radius 


Flexural rigidity 

Shape of undeformed seal ring relative to 
shaft radius 


Viscosity 
Inlet slope 

Lengths shown in Figure A-4 
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Figure A-l is a schematic of a pumping Leningrader seal cross section mounted on 
a shaft. The seal, which is preloaded against the shaft with an interference 
fit, is loaded by a backup spring and by high-pressure gas, which acts behind 
the seal and is separated from the oil by the static sealing land. The long, 
chamfered inlet region provides the pumping action in that it forms a gradual 
inlet which tends to trap or to provide a preferred direction for flow when the 
shaft moves toward the cooling supply. (This is referred to as the forward 
stroke.) During the backstroke, the seal (as drawn in Figure A-l) tends to rub 
and wipe oil away. Thus, any oil that leaks toward the high-pressure gas tends 
to be pumped back toward the cooling oil reservoir. 

Figure A-2 shows the actual dimensions of the seal used in an automotive Stir- 
ling engine. The dotted line denotes the geometry used in the initial modeling 
of the seal discussed below. 

A preliminary model has been developed for a seal with a constant thickness as 
shown in Figure A-3. The pressure profile drawn under the seal is the pressure 
anticipated to occur between the seal and the shaft during the forward stroke. 
The pressure on the far left corresponds to the high-pressure gas, and the pres- 
sure on the right represents the oil pressure. The constant pressure, p^, 
denotes the compressive radial stress associated with the interference fit of 
the ring on the shaft. The gas pressure is assumed to act along the outside of 
the seal with an imaginary secondary seal separating the gas from the oil. The 
deformation of the seal is assumed to be governed by the elasticity equations 
for a thin axisymmetric cylindrical shell. 

D(d 4 w/dx 4 ) + (Et/R 2 )w = -(p - p ) (A-l) 

8 

The geometric variables and coordinate system are shown in Figure A-4; w denotes 
the inward radial deflection, E denotes the elastic modulus of the seal, and D 
denotes the flexural rigidity defined as: 

D = Et^/ [12(1 - v 2 )] 
where v is Poissons ratio for the seal. 
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Fig. A-2 Seal Geometry for Automotive Stirling Engine Application 
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Fig. A-3 Pressures Acting on Pumping Leningrader Seal Used 
in Preliminary Analysis 
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Fig. A-4 Schematic of Geometry Used in the Analysis 
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The quantity, p, enotes the pressure in the interface between the shaft and the 
seal which will be determined by the Reynolds Equation for axisymmetric 
steady-state motion during the forward stroke given by: 

dp/dx = 6yU Q [(h - h*)/h 3 ] (A-2) 

where h is the local hydrodynamic film thickness, y is the viscosity of the oil, 
U q is the average velocity during the forward stroke, and h* is a constant of 
integration to be determined from the boundary conditions. 

The film thickness h is related to the elastic deflection w by: 

h = f(x) - w (A-3) 


where f(x) is the radius of the undeformed seal relative to the shaft radius. 
For the model depicted in Figure A-4, f(x) is given by: 

f -6 , 0 < x < L 

f(x) = | (A-4) 

[_•& - 0x, -L^ < x < 0 


The boundary conditions on the elasticity equation come from the requirement 
that both ends of the seal are free (forces and moments are zero at each end of 
the seal). The boundary conditions on the Reynolds equation comes from the 
requirements that the hydrodynamic pressures equal the prescribed gas pressure 
at the high-pressure side and the prescribed oil pressure on the low-pressure 
side. These may be written as: 

•> v; . 

d^w/dx^ = d 3 w/dx 3 = 0 at x = -L^ and x = L (A-5) 

p = Pg at x = -L^ (A-6) 


p = p at x = L 
r o 


(A-7 ) 


Equations (A-l, Through (A-7) represent a fifth-order nonlinear system corre- 
sponding to the mo ’el under consideration. 
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A number of approximations that are characteristic of elastohydrodynamic lubri- 
cation may be introduced. These are based largely on the fact that the hydrody- 
namic film thickness, h, will be small compared with the interference fit, 6. 
The elastohydrodynamic behavior may be characterized by three zones: 

• an inlet zone where the film develops , 

• a contact zone where the film thickness is nearly constant 

• an exit zone where the pressure drops to ambient. 

This type of behavior has been seen in many other types of elastohydrodynamic 
lubrication, e.g., Hertzian contacts in rolling element bearings, elastomeric 
bearings and foil bearings. 

A.l Contact Zone Solution 


The contact zone is the region extending from x = 0 to the start of the exit zone 
which will be localized near x = L. In that region, from Equation (A-4) , f(x) = 
-6 and from Equation (A-3), w = - 6 - h, which may be approximated by w = 6 since 
| h/ 6 1 « 1. With w being constant, the left-hand terra in Equation (A-l) will 
vanish; hence the relationship 


p - P g = P ± = (Et6)/R 


contact zone pressure 


(A-8) 


Since p as determined in Equation (A-8) is constant in the contact zone, the 
left-hand hand term in Equation (A-2) will vanish and the result will be h = h*, 
constant (contact zone film thickness). 

Equation (A-8) is essentially the dry contact pressure resulting from the inter- 
ference fit. The dry contact profile would also contain a radial shear load at 
the sharp corner at x = 0 . With a small amount of wear, however, this corner 
will be rounded, thus alleviating the radial shear load without radically 
affecting the inlet film shape. Equation (A-8) will then be assumed to prevail 
throughout the contact zone with the value of h* to be determined by matching 
the pressure obtained from the solution to the equations for the inlet zone. 
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A. 2 Inlet Zone Solution 


In order to match the solution in the inlet zone with the contact zone solution, 
the pressure and flow must be continuous at x = 0. This results in the 
constraint 

P = Pg + at x = 0 (A-9) 

and 

dp/dx = 0 at x = 0 
which is equivalent to 

h = h* at x = 0 (A- 10) 


For values of 0 of the order of a few degrees or more, the pressure will vary 
from its limiting values of p^ + p, at x = 0 to p^ over a relatively small 
portion of the inlet zone as sketched in Figure A-3. If the shape of the inlet 
zone is assumed to be unaffected by the pressure distribution in the region, the 
film thickness relationship is 

h = h* - 0x , x < 0 

The stretched variable ? = - 8x/h* may now be introduced into Equation (A-2) to 
obtain 


dp/d? = -.6uU o e/[0h*(l + 5) 3 ] 

with the constraints 

p = p + p . at ? = 0 
g i 

and p = p^ at ? = L^0/h* 


CA-11) 


( A - 1 2 ) 
(A- 13) 
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The quantity L^0/h* will be a large number for practical cases. Equation ( A - 1 3 ) 
may be replaced by its limiting form 

p = p^ at £ ■* •» (A-14) 

Equation (A-ll) may readily be integrated to yield 


6 mU d « r 

* 1ST J 0 oTIT 3 d5 ■ /C9h * ) 

Thus the film thickness h* is given by * 

h* = 3yU /(0p.) = 3viU Q R 2 /(Et68) (A-15) 

which may, in turn, be used to calculate flow. 

A.3 Exit Zone Solution 

The exit zone solution is the region in the neighborhood of x = L where the pres- 
sure must drop from p^ + p^ down to the ambient value of p^. This will, again, 
occur over a relatively short region. The film thickness deflection relation- 
ship for the exit zone as determined from Equations (A-3) and (A-4) is the same 
as that for the contact zone and may be written as 

w = -(h + 6) 

The above relationship may be substituted for w in Equation (A-l). to obtain 
D(d^h/dx^) + Et(5 + h)/R^ = p - p^ 
where x^ = L - x. 

Neglecting h compared with 6 in the second term of the above equation and 
employing Equation (A-8), one obtains 
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(A-16) 


DCdSi/dx^) = p - (p g + P i ) 

One may now differentiate Equation (A-16) with respect to x^ and substitute 
Equation (A-2) for -dp/dx^ to obtain the exit zone equation 

D(d S h/d Xl 5 ) = -6yU (h - h*)/h 3 (A-17) 

1 o 

where h* is now determined from Equation (A-15). 

The boundary conditions at x^ = 0 may be obtained from Equations (A-5), (A-7) 
and (A-16) 

d 2 h/dx x 2 = d 3 h/dx 1 3 = 0 at x 1 = 0 (A-18) 

and 

D(d 4 h/dx 1 4 ) = P Q ‘ P g ' P t at x = 0 (A- 19) 

Without going through the formality of introducing a stretching transformation 
as was done for the inlet zone, the exit zone solution must merge with the 
contact zone solution as x^ **■ <*>. Thus, we look for a bounded solution with h 
h* as Xj 4 «. 

The above system of equations may be solved numerically by Runge-Kutta inte- 
gration, however; if the variation in h is not excessive (|h = h*| « h*) , then 

3 

one may replace h appearing in the denominator of the right-hand side of 

3 

Equation (A-17) with h* to obtain a linear system of equations with constant 
coefficients that may readily be solved algebraicly. 

The algebraic solution is in the form: 

h/h* = 1 + A^e 4 + e ^1^ (A^ cos + A^ sin 

where 

C = [6yU Q /(h* 3 D)3 1/5 x x , X = cos (2 tt/ 5) and X 2 = sin (2ir/5). 
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The constants A^, and A^ are readily determined from Equations (A-18) and 
(A-19) . 

A.4 Sample Solution 

Film thickness and pressure profiles have been calculated corresponding to the 
geometry shown in Figure A-2 with the trapezoidal section modeled by the dotted 
line shown in the figure. A fluid viscosity of 55 cp, an elastic modulus of 
1.72 GPa, and a Poisson's ratio of 0.41 have been used in the computation. 

As shown on the scale in Figure A-5 , the inlet zone film thickness profile is 
very steep, indicating that only the portion of the inlet zone very near the 
contact zone is relevant in affecting the film thickness profile, and that the 
major portion of the seal to the left of this region is probably not necessary. 
This is again shown in Figure A-6, where the pressure, approximately 1 mm before 
the start of the contact zone, is very nearly that of the gas (10 MPa). The 
pressure rises very steeply and joins the pressure in the contact zone smoothly, 
although it appears to have a sharp corner on the scale in Figure A-6. The pres- 
sure profile in the contact zone is constant and equal to the sum of the gas 
pressure plus the interference pressure (the pressure that would be associated 
with interference fit statically). Near the exit region, the pressure starts to 
oscillate, peaks to 18 MPa, and then falls rapidly to 0 (the oil pressure). 

Examination of the rapid falloff in pressure at the exit region shows that if 
the exit oil pressure were raised to be equal to, or even greater than, the gas 
pressure, there would be little change in the pressure profile except in the 
immediate vicinity of the exit zone. The film thickness in the contact zone 
would be virtually unaffected. This indicates that the ambient pressure differ- 
ence has little effect on the film thickness or the flow rate since the flow is 
completely dominated by viscous forces. Thus, the seal could be made to provide 
a film thickness in either direction by adding a chamfer on the reverse side. 
Such chamfers, which were recommended and tested as part of the Automotive Stir- 
ling Engine Program, have been shown to be effective in providing lubrication 
without excessive leakage in either direction. 
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FILM THICKNESS (MICRONS) 



PRESSURE (MPA) 






APPENDIX B 


COMPUTER PROGRAM "RING" 

The computer program RING which has been used throughout the report is listed 
herein. It is written in IBM FORTRAN 77 (VS), and uses an IMSL subroutine 
ZSCNT. A FORTRAN IV Listing of ZSCNT and its subroutines is also included for 
completeness. The FORTRAN listings are preceded with an input description, an 
output description, and sample inputs and outputs. 
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B.l INPUT DESCRIPTION 


AIL input variables are on a namelist (NAMELIST/INPUTS/) . Symbols in the left 
coLumn denote the FORTRAN names of the namelist variables. Symbols appearing 
on the right correspond to the nomenclature given at the beginning of the 
report. 

Input Variable Definitions 

STRING A character string of up to 60 characters to identify job 
IDIM Type of input 

0 - Dimensionless (default value) 

1 - Dimensional (any consistent set of units) 

Dimensional Inputs (used only if IDIM = 1) 

NI Frequency (Hz), f 

Si Stroke (L), s 

Cl Clearance (L), c 

ELI Bearing length (L), L 

EL1I Length of nonbearing portion of ring (L), Li 

El Length over which p 0 acts (L), e 

RI Shaft radius (L), R 

TI Ring thickness (L), t 

ClI Initial taper, 5 (defaults to 0) 

POI Preload pressure (F/L^), p 0 

PFI Reservoir pressure (F/L^), pf 

RMUI Lubricant viscosity (FT/L^), U 

EMOD Ring modulus of elasticity (F/L^), E 

POIS Ring Poisson's Ratio, v 

Nondimensional Inputs 


a 

B 


AL 

BET 

EPS 


e 
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ELI 

Cl 

PO 

PF 

ELOS 

NDX 

NDP 

ICAV 

I PR 

NMAX 

NSIG 

XF 


1*1 

Sl/C (defaults to 0) 

Po 

Pf 

L/s for starvation calculation (defaults to 0) 

A vector of length 3 containing number of increments for 
Runge-Kutta integration 

NDX(l) = Number of increments inl-e<5<l (defaults to 40) 
NDX(2) = Number of increments inQ<5<l“ e (defaults to 25) 
NDX(3) = Number of increments in -Li < C 5 ® (defaults to 25) 

A vector of length (3) containing number of points at which 
pressures etc. are saved for printout. NDX(i) should be inte- 
gral multiples of NDP(i). (Defaults are NDP(l) = 10, NDP(2) = 
25, NDP(3) = 25.) 

Flag calculation of cavitation (backstroke) 

0 - No cavitation (negative pressures allowed) 

2 - Cavitation included, no negative pressures (default value) 

Print Flag 

0 - Short output (default value) 

1 - Longer output (print pressure profile etc.) 

Maximum iterations for secant method solution (default is 10) 

Number of significant digits in accuracy of solution using 
secant method (default is 4) 

A vector of length (2) containing initial guesses for forward 

stroke secant solution 

XF(1) = Initial guess for h at ( = 1 

XF(2) = Initial guess for dh/d£ at £ = 1 
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If XF(1) and XF(2) are both set equal to 0, the program will 
attempt to find its own initial guesses. In subsequent runs, 
values determined from the previous run are used unless 
explicitly specified. 

A vector length (2) containing initial guesses for h and dh/dij 
for reverse stroke secant solution. Same rules apply to XB as 


B. 2 OUTPUT DESCRIPTION 


The output starts with a listing of the inputs which are defined in D.l. 
Dimensional inputs are listed only if IDIM = 1. All outputs are dimension- 
less. 

The Elastic Solution refers to the deflection of the ring when hydrodynamic 
forces are not present. If the ring is unclampied, h(l) > 0, one line of 
output is given for the Elastic Solution. If the ring is clamped, h(l) < 0, 
two lines of output are given. The first line refers to the ring deflection 
relative to an imaginary shaft that would occur if the shaft were not present 
(the negative value of h noting the interference). The second line of output 
includes the effect of the interference force exerted by the shaft on the ring 
at ^ = 1 when h(l) = 0. 


Thw Pumping Flow Solution and Back Flow Solution refer to forward and reverse 
strokes respectively. They will contain one line of output each at IPR = 0 or 
multiple outputs if IPR = 1. 


The following table relates the output caption with the symbols for variables 
as defined in the nomenclature for the elastic, pumping flow and back flow 
solution. 


Output Caption 


Algebraic Symbol 


X 

H 

H' 

H" 

H'" 

PRES 

FORCE 


K 

h 

dh/d? 

d 2 h/d^ 2 

d 3 h/d£ 3 

P 

F 


The final output listed under flows relates to dimensionless flow and cavita- 
tion output and are given as follows: 
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Output Caption 


Algebraic Symbol 


PUMPING 

BACK 

XCAV 

NET FLOW 
STARVED NET 


FLOW 


Kp (Forward Stroke) 
Kr (Reverse Stroke) 

Kp - Kr 
Kf eff ' Kr 


l?.l 


B. 3 SAMPLE INPUT AND OUTPUT 


The sample input and output contained herein correspond to three solutions for 
a babbitt pumping ring with loading pressures, p 0 , of 3.45, 5.17 and 6.89 MPa 
(500, 750 and 1000 psi). A full printout (IPR = 1) was requested at p 0 = 3.45 
MPa (500 psi). It should be noted that the ring is predicted to be clamped 
during the backstroke at p 0 = 6.89 MPa (1000 psi). 
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FILE: RING INPUT A 1 MTI FRI 06/17/83 13:31:29 


& INPUTS 

STRING=' BABBIT RING A-l-A-1 P0=500 COMPLETE OUTPUT (IPR=1) ' 

IDIM=1 , NI =35., SI = 2., C1I =0., 

Cl =.5E-3, RMUI =. 885E-5 , ELI =.267, . POI =500., 

PFI =0., El =.115, ELI I =.298, RI =.375, 

TI =.047, EMOD =7.5E6, ' POIS =.36, 

XF = 0. ,0. , XB = 0 . , 0 . , IPR=1 , 

&END 

&INPUTS P0I=750. , STRING =' BABBIT RING A-l-A-l P0=750' ,IPR=0,&END 
&INPUTS PO 1=1000. .STRING =' BABBIT RING A-l-A-1 P0=1000 ' ,&END 
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* PUMPING RING ANALYSIS PROGRAM * 


INPUTS 


icictcicicicie'kic'k'kit'kifieicie 


STRING 

! = BABBIT RING A-l 

-A-l 

P0=500 COMPLETE 01 

DIMENSIONAL INPUTS 



ELI 

= 0.26700000D+00; 

EL1I 

= 0.29800000D+00 

El 

= 0.11500000D+00; 

Cl 

= 0.50000000D-03 

C1I 

= 0.00000000D+00; 

TI 

= 0.47000000D-01 

SI 

= 0.20000000D+01; 

RI 

= 0 . 37500000D+00 

RMUI 

= 0.88500000D-05; 

NI 

= 0.35000000D+02 

POI 

= 0.50000000D+03; 

PFI 

= O.OOOOOOOOD+OO 

EMOD 

» 0. 75000000D+07 ; 

POIS 

= 0.36000000D+00 

<SrHnWnlk 

IDIM 

= i; 


'ricirick-irk'k-ifk'kit'ificicick-i''. 

AL 

= 0.66085654D-02; 

BET 

= 0.71535454D+01; 

EPS 

= 0.43071161D+00; 

PO 

= 0.62976163D-01; 

ELI 

= 0. 11161049D+01 ; 

PF 

= O.OOOOOOOOD+OO; 

STARV 

= 0.32041966D+00; 

Cl 

= O.OOOOOOOOD+OO; 




NDX 

NDP 

NMAX 

IPR 

XF 

XB 


40, 

10 , 


25, 

25, 


25; 

25; 


10; NSIG = 

1; icav = 

O.OOOOOD+OO 

0.00000D+00 


4; 

2 ; 

O.OOOOOD+OO 

O.OOOOOD+OO 


OUTPUTS 


* ELASTIC SOLUTION * 


X 

,0000 

H 

0.48845 -0. 

H' H' ’ 

67281 0.00000 

H t " 

0.00000 

PRES 

0.00000 

FORCE 

0.00000 



* PUMPING FLOW 

SOLUTION * 



X 

H 

H' 

H" 

H" ' 

PRES 

1.0000 

0.81814 

-0.69559 

0.00000 

0.00000 

0.00000 

0.9914 

0.82413 

-0.69558 

-0.00147 

0.33645 

0.00375 

0.9742 

0.83611 

-0.69548 

-0.01254 

0.93187 

0.01074 

0.9569 

0.84809 

-0.69510 

-0.03304 

1.43418 

0.01709 

0.9397 

0.86006 

-0.69429 

-0.06149 

1.85487 

0.02285 

0.9225 

0.87201 

-0.69294 

-0.09655 

2.20463 

0.02805 

0.9052 

0.88393 

-0.69093 

-0.13710 

2.49334 

0.03274 

0.8880 

0.89582 

-0.68819 

-0.18216 

2.73019 

0.03695 

0.8708 

0.90764 

-0.68464 

-0.23092 

2.92369 

0.04071 

0.8536 

0.91940 

-0.68022 

-0.28270 

3.08173 

0.04405 

0.8363 

0.93108 

-0.67488 

-0.33695 

3.21161 

0.04701 

0.8191 

0.94265 

-0.66859 

-0.39324 

3.32008 

0.04960 

0.8019 

0.95411 

-0.66132 

-0.45126 

3.41339 

0.05186 
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0.7846 

0.96543 

-0.65304 

-0.51080 

3.49729 

0.05380 

0.7674 

0.97660 

-0.64371 

-0.57175 

3.57709 

0.05545 

0.7502 

0.98761 

-0.63333 

-0.63406 

3.65768 

0.05682 

0.7330 

0.99842 

-0.62186 

-0.69781 

3.74354 

0.05793 

0.7157 

1.00903 

-0.60927 

-0.76311 

3.83880 

0.05881 

0.6985 

1.01941 

-0.59555 

-0.83016 

3.94719 

0.05946 

0.6813 

1.02954 

-0.58066 

-0 .89921 

4.07216 

0.05990 

0.6640 

1.03941 

-0.56456 

-0.97058 

4.21679 

0.06015 

0.6468 

1.04899 

-0.54720 

-1.04464 

4.38389 

0.06021 

0.6296 

1.05825 

-0.52854 

-1.12178 

4.57598 

0.06011 

0.6124 

1.06719 

-0.50853 

-1.20247 

4.79530 

0.05984 

0.5951 

1.07577 

-0.48709 

-1.28718 

5.04382 

0.05943 

0.5779 

1.08396 

-0.46415 

-1.37644 

5.32328 

0.05887 

0.5693 

1.08791 

-0.45209 

-1.42295 

5.47508 

0.05855 

0.5579 

1.09297 

-0.43555 

-1.48207 

4.91228 

0.05807 

0.5351 

1.10249 

-0.40062 

-1.58150 

3.83132 

0.05695 

0.5124 

1.11120 

-0.36370 

-1.65702 

2.81167 

0.05564 

0.4896 

1.11904 

-0.32532 

-1.71003 

1.85497 

0.05416 

0.4668 

1.12600 

-0.28598 

-1.74199 

0.96235 

0.05252 

0.4440 

1.13206 

-0.24614 

-1.75435 

0.13448 

0.05073 

0.4213 

1.13721 

-0.20622 

-1.74861 

-0.62831 

0.04881 

0.3985 

1.14146 

-0.16663 

-1.72623 

-1.32603 

0.04676 

0.3757 

1.14481 

-0.12772 

-1.68871 

-1.95894 

0.04461 

0.3530 

1.14728 

-0.08982 

-1.63751 

-2.52748 

0.04236 

0.3302 

1.14891 

-0.05323 

-1.57408 

-3.03230 

0.04001 

0.3074 

1.14972 

-0.01821 

-1.49988 

-3.47412 

0.03759 

0.2846 

1.14975 

0.01500 

-1.41633 

-3.85377 

0.03509 

0.2619 

1.14905 

0.04623 

-1.32484 

-4.17212 

0.03252 

0.2391 

1.14766 

0.07529 

-1.22678 

-4.43005 

0.02989 

0.2163 

1.14564 

0.10206 

-1.12353 

-4.62844 

0.02721 

0.1936 

1.14303 

0.12643 

-1.01643 

-4.76813 

0.02448 

0.1708 

1.13990 

0.14833 

-0.90681 

-4.84992 

0.02171 

0.1480 

1.13630 

0.16772 

-0.79599 

-4.87453 

0.01890 

0.1252 

1.13228 

0.18459 

-0.68524 

-4.84260 

0.01606 

0.1025 

1.12791 

0.19894 

-0.57586 

-4.75471 

0.01319 

0.0797 

1.12324 

0.21083 

-0.46912 

-4.61129 

0.01029 

0.0569 

1.11833 

0.22034 

-0.36627 

-4.41272 

0.00737 

0.0342 

1.11322 

0.22755 

-0.26857 

-4.15924 

0.00443 

0.0114 

1.10798 

0.23262 

-0.17726 

-3.85102 

0.00148 

0.0000 

1.10532 

0.23439 

-0.13439 

-3.67640 

0.00000 



* BACK FLOW 

SOLUTION * 




X 

H 

H' 

H" 

H”’ 

PRES 

1.0000 

0.48845 

-0.67281 

0.00000 

0.00000 

0.00000 

0.9914 

0.49424 

-0.67281 

0.00033 

-0.07580 

0.00000 

0.9742 

0.50584 

-0.67284 

0.00279 

-0.20473 

0.00000 

0.9569 

0.51743 

-0.67292 

0.00721 

-0.30344 

0.00000 

0.9397 

0.52902 

-0.67309 

0.01307 

-0.37193 

0.00000 

0.9225 

0.54062 

-0.67338 

0.01985 

-0.41019 

0.00000 

0.9052 

0.55223 

-0.67378 

0.02703 

-0.41820 

0.00000 

0.8880 

0.56384 

-0.67431 

0.03409 

-0.39595 

0.00000 

0.8708 

0.57546 

-0.67495 

0.04050 

-0.34341 

0.00000 

0.8536 

0.58710 

-0.67570 

0.04575 

-0.26055 

0.00000 

0.8363 

0.59875 

-0.67652 

0.04931 

-0.14735 

0.00000 
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0.8191 

0.61041 

-0.67738 

0.05065 

-0.00376 

0.00000 

0.8019 

0.62209 

-0.67825 

0.04926 

0.17025 

0.00000 

0.7846 

0.63378 

-0.67906 

0.04461 

0.37472 

0.00000 

0.7674 

0.64548 

-0.67976 

0.03617 

0.60970 

0.00000 

0.7502 

0.65720 

-0.68028 

0.02343 

0.87520 

0.00000 

0.7330 

0.66892 

-0.68054 

0.00584 

1.17126 

0.00000 

0.7157 

0.68065 

-0.68045 

-0.01711 

1.49788 

0.00000 

0.6985 

0.69237 

-0.67992 

-0.04595 

1.85507 

0.00000 

0.6813 

0.70407 

-0.67883 

-0.08120 

2.24279 

0.00000 

0.6640 

0.71575 

-0.67708 

-0.12340 

2.66099 

0.00000 

0.6468 

0.72740 

-0.67454 

-0.17307 

3.10961 

0.00000 

0.6296 

0.73899 

-0.67107 

-0.23072 

3.58851 

0.00000 

0.6124 

0.75051 

-0.66654 

-0.29689 

4.09756 

0.00000 

0.5951 

0.76195 

-0.66079 

-0.37208 

4.63653 

0.00000 

0.5779 

0.77327 

-0.65366 

-0.45682 

5.20518 

0.00000 

0.5693 

0.77889 

-0.64953 

-0.50292 

5.50054 

0.00000 

0.5579 

0.78625 

-0.64346 

-0.56341 

5.12594 

0.00000 

0.5351 

0.80075 

-0.62936 

-0.67194 

4.41448 

0.00000 

0.5124 

0.81490 

-0.61297 

' -0.76483 

3.75238 

0.00000 

0.4896 

0.82865 

-0.59464 

-0.84320 

3.13837 

0.00000 

0.4668 

0.84196 

-0.57467 

-0.90812 

2.57100 

0.00000 

0.4440 

0.85481 

-0.55337 

-0.96063 

2.04872 

0.00000 

0.4213 

0.86716 

-0.53101 

-1.00175 

1.56985 

0.00000 

0.3985 

0.87899 

-0.50783 

-1.03244 

1.13265 

0.00000 

0.3757 

0.89028 

-0.48406 

-1.05364 

0.73528 

0.00000 

0.3530 

0.90103 

-0.45991 

-1.06622 

0.37590 

0.00000 

0.3302 

0.91123 

-0.43556 

-1.07103 

0.05260 

0.00000 

0.3074 

0.92087 

-0.41118 

-1.06887 

-0.23653 

0.00000 

0.2846 

0.92995 

-0.38693 

-1.06050 

-0.49338 

0.00000 

0.2619 

0.93849 

-0.36293 

-1.04663 

-0.71988 

0.00000 

0.2391 

0.94649 

-0.33930 

-1.02793 

-0.91790 

0.00000 

0.2163 

0.95395 

-0.31614 

-1.00503 

-1.08929 

0.00000 

0.1936 

0.96089 

-0.29355 

-0.97851 

-1.23587 

0.00000 

0.1708 

0.96732 

-0.27160 

-0.94892 

-1.35941 

0.00000 

0.1480 

0.97326 

-0.25036 

-0.91676 

-1.46164 

0.00000 

0.1252 

0.97873 

-0.22987 

-0.88250 

-1.54421 

0.00000 

0.1025 

0.98374 

-0.21018 

-0.84657 

-1.60875 

0.00000 

0.0797 

0.98831 

-0.19132 

-0.80936 

-1.65679 

0.00000 

0.0569 

0.99246 

-0.17332 

-0.77123 

-1.68981 

0.00000 

0.0342 

0.99621 

-0.15620 

-0.73251 

-1.70922 

0.00000 

0.0114 

0.99958 

-0.13996 

-0.69348 

-1.71637 

0.00000 

-0.0000 

1.00113 

-0.13218 

-0.67394 

-1.71575 

0.00000 


* FLOWS * 

PUMPING BACK XCAV NET FLOW STARVED NET FLOW 

0 . 106219D+01 0.488604D+00 0.999770D+00 0.573582D+00 0.296125D+00 

********'ir^************************************Vf****************'**-****** 


126 



* Pl/MPING RING ANALYSIS PROGRAM * 


INPUTS 


ieirirtrln 


kkkkirkklrkkirkkkkkkkkirlrtrtrirfrtr 

■kirklrklriri 

i.. 1 . j.j..* -< » ■ > « 

IXWn n KHJf 

STRING = BABBIT RING A-l-A-1 

P0=750 

i-irfrk-k-k-k-k-k-k-k-k-k-k-k-klc-k-k-k-l.-k-k-k-k-k 

i . i f . - > i 

HfiSr******* 

DIMENSIONAL INPUTS 




ELI 

El 

CII 

SI 

RMUI 

= 0.26700000D+00; EL1I 
= 0.11500000D+00; Cl 
= O.OOOOOOOOD+OO; TI 
= 0.20000000D+01; RI 
= 0.88500000D-05; NI 

= 0.29800000D+00; 
= 0.50000000D-03; 
= 0.47000000D-01; 
= 0.37500000D+00; 
= 0.35000000D+02; 



POI 

EMOD 

= 0. 75000000D+03; PFI 
= 0. 75000000D+07 ; POIS 

= 0.00000000D+00; 
= 0.36000000D+00; 

^ ^ ^ ^ ^ y 


IDIM 

AL 

EPS 

ELI 

STARV 

= i; 

= 0.66085654D-02; BET 
= 0.43071161D+00; PO 
= 0.11161049D+01; PF 
= 0.32041966D+00; Cl 

= 0 . 71535454D+01 ; 
= 0 . 94464244D-01 ; 
= O.OOOOOOOOD+OO; 
= O.OOOOOOOOD+OO; 



NDX 

NDP 

= 40, 25, 25; 

= 10, 25, 25; 




NMAX 
I PR 
XF 
XB 

•kirkirki 

= 10; NSIG = 4; 

= 0; ICAV = 2; 

= 0.81814D+00 -0.69559D+00 

= 0.48845D+00 -0.67281D+00 

trfrfc1cltirir£ 

rtnWrtrWnHrtf 


OUTPUTS 



* ELASTIC SOLUTION * 

X H H' H" H"’ PRES FORCE 

1.0000 0.23267 -1.00922 0.00000 0.00000 0.00000 0.00000 


* PUMPING FLOW SOLUTION * 


X 

H 

H' 

H' ' 

H" ' 

PRES 

1.0000 

0.71873 

-0.93560 

0.00000 

0.00000 

0.00000 


* BACK FLOW SOLUTION * 


X 

H 

H' 

H' ' 

H’” 

PRES 

1.0000 

0.23267 

-1.00922 

0.00000 

0.00000 

0.00000 


* FLOWS * 


PUMPING BACK XCAV NET FLOW STARVED NET FLOW 

0.100210D+01 0.232747D+00 0.999927D+00 0.769350D+00 0.329915D+00 
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* PUMPING RING ANALYSIS PROGRAM * 


INPUTS 

********************** *** ****************** *** ****-********************** 


STRING = BABBIT RING A-l-A-1 P0=1000 


************** 


rtf*********-. 




rtrtf ************* * * ************ 


DIMENSIONAL INPUTS 


ELI = 0.26700000D+00 
El = 0 . 1 1500000D+00 
C1I = O.OOOOOOOOD+OO 
SI = 0 . 20000000D+01 
RMUI = 0.88500000D-05 


EL1I = 0.29800000D+00; 
Cl = 0.50000000D-03; 
TI = 0.47000000D-01; 
RI = 0 . 37500000D+00 ; 
NI = 0.35000000D+02; 


POI = 0 . 10000000D+04; PFI = O.OOOOOOOOD+OO; 

EMOD = 0.75000000D+07; POIS = 0.36000000D+00; 


t ************ 


IDIM ■ 

is 





AL 

0.66085654D-02 

; BET 

= 0.71535454D+01; 



EPS 

0.43071161D+00 

; PO 

= 0.12595233D+00; 



ELI 

0. 11161049D+01 

; PF 

= O.OOOOOOOOD+OO; 



STARV = 

0.32041966D+00 

; Cl 

= O.OOOOOOOOD+OO; 



NDX 

40, 25, 

25; 




NDP 

10, 25, 

25; 




NMAX = 

10; NSIG = 

4; 




IPR 

0; ICAV = 

2; 




XF 

0.7 1873D+00 

-0 

.93560D+00 



XB 

ie'ieirtricieic'i 

0.23267D+00 

-0 

*****- 

. 10092D+01 

<■-************ ********* 

********** 

rtf************ 

OUTPUTS 


'k'kic'lcJc, 

********************** 

********** 

c'k'tc'k'ic-k'k'k'k'k'k'k'k'k 



* 

ELASTIC SOLUTION * 



X 

H 

H’ 

H" H"’ 

PRES 

FORCE 

1.0000 

-0.02310 -1. 

34562 

0.00000 0.00000 

0.00000 

0.00000 

1.0000 

0.00000 -1. 

28833 

0.00000 -0.70479 

0.00000 

0.00065 



* PUMPING FLOW SOLUTION * 



X 

< H 


H' H" 

H'" 

PRES 

1.0000 0.63082 

-1 

.13187 0.00000 

0.00000 

0.00000 



* ] 

BACK FLOW SOLUTION * 



X 

H 

H' 

H” H"' 

PRES 

FORCE 

1.0000 

0.00000 -1. 

28833 

0.00000 -0.70479 

0.00000 

0.00065 


* FLOWS * 


PUMPING BACK XCAV NET FLOW STARVED NET FLOW 

0.929129D+00 0.000000D+00 0.100000D+0I 0.929129D+00 0.333707D+00 


********************************* 
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B. 4 FORTRAN LISTING 


FILE: RING FORTRAN L4 MTI 


MON 01/20/86 14:39:26 


PAGE 1 OF 38 


PROGRAM RING 


XLAM=.739*ELOS**-.585 IS NOW BUILT INTO PROGRAM 


C$ VS 

C> 

CNOTEJ 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FUNCTION 


RESTRICTIONS 


REMARKS 


- PUMPING RING ANALYSIS PROGRAM 

- STEADY- STATE ANALYSIS ONLY 

- VS FORTRAN (FORTRAN 77) 

IMPLICIT DOUBLE PRECISION REAL*8 


EXTERNAL REFERENCES 
FORTRAN ROUTINES 

IMSL ROUTINES 


USER ROUTINES 


DATAN 

ZSCNT ; SOLVES THE SYSTEM OF NON-LINEAR 
EQUATIONS ; LISTING PROVIDED HEREIN BY 
PERMISSION OF IMSL. 


AIN ; 
AMU ; 
CALCD 

CHECK 

CONST 

CONST2 
DFN1 
DFN2 
DFN3 
ELAS 

ERRMSG 


COMPUTE INVERSE OF 2X2 MATRIX 
MULTIPLY 2 2X2 MATRICES 
; CALCULATE ELASTIC INFLUENCE 
COEFFICIENTS C,D 
; FOR MULTIPLE RUNS CHECK INPUTS 
FOR RECALCULATION OF C,D 
; CALCULATE SLIDER BEARING PRESSURE 
CONSTANTS 
; SAME AS ABOVE 

DERIVATIVE FUNCTION USED BY RUK 


DETERMINE ELASTIC SOLUTION 
(NO HYDRODYNAMICS) 

; PRINT ERROR MESSAGE IF *ZSCNT* 
CONVERGED 


RIN00010 
RIN00020 
-RIN00030 
RIN00040 
RIN00050 
RIN00060 
RIN00070 
RIN00080 
RIN00090 
RIN00100 
RIN00110 
RIN00120 
RIN00130 
RIN00140 
RIN00150 
RIN00160 
RIN00170 
RIN00180 
RIN00190 
RIN00200 
RIN00210 
RIN00220 
RIN00230 
RIN00240 
RIN00250 
RIN00260 
RIN00270 
RIN00280 
RIN00290 
RIN00300 
RIN00310 
RIN00320 
RIN00330 
RIN00340 
RIN00350 
NOTRIN00360 
RIN00370 


EVAL ; DEFINE NON-LINEAR SYSTEM IN H AND H'RIN00380 


INPUT/OUTPUT: 
UNIT 

4 

5 

6 


TO BE SOLVED BY *ZSCNT* 

P ; PRESSURE FUNCTION 

PRT ; CONVERT FROM W TO H AND PRINT 

PRTOUT ; PRINT OUT RESULTS 

RUK : RUNGE-KUTTA 


DESCRIPTION 
TERMINAL I/O 

INPUT FILE IN NAMELIST FORMAT 
OUTPUT FILE 


INPUT VARIABLE DEFINITIONS 

* NOTE : (D) INDICATES VARIABLE 

NAME DESCRIPTION 

NAMELIST /INPUTS/ 


HAS A DEFAULT VALUE 


RIN00390 

RIN00400 

RIN00410 

RIN00420 

RIN00430 

RIN00440 

RIN00450 

RIN00460 

RIN00470 

RIN00480 

RIN00490 

RIN00500 

RIN00510 

RIN00520 

RIN00530 

RIN00540 

RIN00550 
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FILE: RING FORTRAN L4 mti 


MON 01/20/86 14:39:26 


PAGE 2 OF 38 


c 



RIN00560 

c 

STRING 

CHARACTER STRING TO IDENTIFY JOB (MAX. 60 CHARS.) 

RIN00570 

c 

IDIM (D) 

TYPE OF INPUT 

RIN00580 

c 


0 - DIMENSIONLESS 

RIN00590 

c 


1 - DIMENSIONAL (ANY CONSISTENT SET OF UNITS) 

RIN00600 

c 



RIN00610 

c 

DIMENSIONAL INPUTS (NOTE: U0=2.*NI*SI) 

RIN00620 

c 


■ 

RIN00630 

c 

NI 

FREQUENCY (HZ) 

RIN00640 

c 

SI 

STROKE (L) 

RIN00650 

c 

Cl 

CLEARANCE (L) 

RIN00660 

c 

ELI 

BEARING LENGTH (L) 

RINQQ670 

c 

EL1I 

LENGTH OF NON-BEARING PORTION OF RING (L) 

RIN00680 

c 

El 

LENGTH FROM END OF BEARING PRELOADED WITH POI (L) 

RIN00690 

c 

RI 

RING RADIUS (L) 

RIN00700 

c 

TI 

RING THICKNESS (L) 

RIN00710 

c 

C1I (D) 

BEARING SLOPE 

RIN00720 

c 

POI 

PRELOAD PRESSURE (F/L**2) 

RIN00730 

c 

PFI 

RESERVOIR PRESSURE (F/L**2) 

RIN00740 

c 

RMUI 

LUBRICANT VISCOSITY (FT/L**2) 

RIN00750 

c 

EMOD 

RING MODULUS OF ELASTICITY (F/L**2) 

RIN00760 

c 

POIS 

RING POISSON'S RATIO 

RIN00770 

c 



RIN00780 

c 

NON-DIMENSIONAL INPUTS 

RIN00790 

c 



RIN00800 

c 

AL 

(TI*RI )**2/( 12*ELI**4*( 1 .-POIS)**2) ) 

RIN00810 

c 

BET 

(T*RMUI*UO*RI**2*ELI)/(CI**3*TI*EMOD) 

RIN00820 

c 

EPS 

LENGTH FROM END OF PRELOAD PO (El /ELI) 

RIN00830 

c 

ELI 

LENGTH OF NON-BEARING PORTION OF RING (ELll/ELI) 

RIN00840 

c 

ELOS (D) 

L/S LAND TO STROKE RATIO FOR STARVATION CALCULATIONRIN00850 

c 

XLAM (D) 

MULTIPLIES ELOS FOR INCREASED STARVATION 

RIN00860 

c 

Cl (D) 

SLOPE OF BEARING 

RIN00870 

c 

PO 

PRELOAD PRESSURE (CI**2/ (6*RMUI*U0*ELI ) )*POI 

RIN00880 

c 

PF 

RESERVOIR PRESSURE (CI**2/(6*RMUI*U0*ELI ) )*PFI 

RIN00890 

c 



RIN00900 

c 

NDX(3) (D) 

DELTA X INCREMENTS FOR RUK 

RIN00910 

c 

NDP(3 ) (D) 

tl H fl It It 

RIN00920 

c 



RIN00930 

c 

ICAV (D) 

FLAG FOR CALCULATION OF CAVITATION (BACK STROKE) 

RIN00940 

c 


0 - NO CAVITATION (NEG. PRESSURES ALLOWED) 

RIN00950 

c 


1 - CAVITATION (NO NEG. PRESSURES) 

RIN00960 

c 


2 - FIND ITS OWN SOLUTION 

RIN00970 

c 

IPR (D) 

PRINT FLAG 

RIN00980 

c 


0 - SHORT OUTPUT 

RIN00990 

c 


1 - LONGER OUTPUT (PRINT PRESSURE PROFILE) 

RIN01000 

c 

NMAX (D) 

MAX. ITERATIONS FOR *ZSCNT* 

RIN01010 

c 

NSIG (D) 

NO. OF SIGNIFICANT DIGITS IN ACCURACY 

RIN01020 

c 


OF SOLUTION USING *ZSCNT* 

RIN01030 

c 

XF(2) 

INITIAL GUESS FOR H AND H’ FOR FORWARD STROKE 

RIN01040 

c 

XB(2) 

INITIAL GUESS FOR H AND H' FOR BACKWARD STROKE 

RIN01050 

c 



RIN01060 

c 



RIN01070 




“KlflJUlUoU 


IMPLICIT REAL*8 (A-H.O-Z) 

RIN01090 

c 



RIN01100 
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non o o n n non 


FILE: RING FORTRAN L4 MTI 


MON 01/20/86 14:39:26 


PAGE 3 OF 38 


COMMON/BDIM /Cl , ELI , CONI , IDIM 
COMMON/BPAR /AL,BET,P0,EPS,PF 

COMMON/BCD /XLOC(IOO) ,C( 100 ) ,D(100) ,NDX(3) ,NDP(3) ,DX(3) ,DP(3) 
COMMON/ BCOE FF / RK , RC , R 1 , R2 , XCAV , I CAV 
COMMON/BINT /IELAS,IBACK,IPR 
COMMON/BLAST /X1LAST,X2LAST,F1LAST, F2LAST 
COMMON/BELAS / FORCE , HELAS , HI ELAS , H2 ELAS , H3ELAS , W1 ( 4 ) , W2 ( 4 ) 
COMMON/ BPROF /Cl 
COMMON/BFLAG /IFLAG,IPL0AD 

DIMENSION PAR(1),XX(2),XF(2),XB(2),WK(68) 

REAL*8 NFL0W,NI 
CHARACTER*60 STRING 
LOGICAL RECALC 
EXTERNAL EVAL , EVAL3 , EVAL4 
NAMELIST /INPUTS/ 

DIMENSIONAL INPUTS 

♦ STRING, NI , SI ,C1I , Cl ,RMUI , ELI ,P0I ,PFI ,EI ,EL1I ,RI ,TI ,EMOD,POIS 

NON-DIMENSIONAL INPUTS 

+ , BET , PO , EPS , AL , ELI , ELOS , XLAM , 
+NDX,NDP,PF,NMAX,IPR,NSIG,ICAV,XF,XB,C1,IDIM 

PI=4.*DATAN( 1.D0) 

DEFAULTS 

IDIM=0 
ELOS=0. 

XLAM=1. 

C1I=0, 

C1=0. 

NDX( l)=40 
NDX(2)=25 
NDX(3)=25 
NDP( 1)=10 
NDP(2)=25 
NDP(3)=25 
NMAX=10 
NSIG=4 
ICAV=2 
IPR=0 
C 

IRUN=0 

2001 READ(05, INPUTS, END=2000) 

C 

IRUN=IRUN+1 

WRITE(6,'("1* PUMPING RING ANALYSIS PROGRAM *",/)') 
WRITE(6,'(" INPUTS")') 

WRITE(6,'(1X,72("*"))') 

WRITE(6,'(" STRING = ",A60)') STRING 
WRITE(6,'(1X,72("*"))') 


RIN01110 
RIN01120 
RIN01130 
RIN01140 
RIN01150 
RIN01160 
RIN01170 
RIN01180 
RIN01190 
RIN01200 
RIN01210 
RIN01220 
RIN01230 
RIN01240 
RIN01250 
RIN01260 
RIN01270 
RIN01280 
RIN01290 
RIN01300 
RIN01310 
RIN01320 
RIN01330 
RIN01340 
RIN01350 
RIN01360 
RIN01370 
RIN01380 
RIN01390 
RIN01400 
RIN01410 
RIN01420 
RIN01430 
RIN01440 
RIN01450 
RIN01460 
RIN01470 
RIN01480 
RIN01490 
RIN01500 
RIN015 10 
RIN01520 
RIN01530 
RIN01540 
RIN01550 
RIN01560 
RIN01570 
RIN01580 
RIN01590 
RIN01600 
RIN01610 
RIN01620 
RIN01630 
RIN01640 
RIN01650 
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FILE: RING FORTRAN L4 MTI 

MON 01/20/86 14:39:26 PAGE 4 OF 38 

IF(IDIM.NE.O)THEN 



RIN01660 

WRITE(6,'(” DIMENSIONAL INPUTS",/)') 


RIN01670 

WRITE(6,'(" ELI 

= ",E15.8,"; EL1I 

= ' ' , E15 .8, ";"•)') 

RIN01680 

+ELI ,EL1I 



RIN01690 

WRITE(6,'(” El 

= ",E15.8,"; Cl 

= ' ' ,E15.8," ; " ) ' )EI 

,CIRIN01700 

WRITE(6, '(" ClI 

= ",E15.8,"; TI 

= ' ' ,E15.8, ";")') 

RIN01710 

+C1I ,TI 



RIN01720 

WRITE(6,'(" SI 

= ",E15.8,"; RI 

=»' ' ,E15.8, " ; " ) ' )SI 

,RIRIN01730 

WRITE(6,'(" RMUI 

= ' ' ,E15.8," ; NI 

= ",E15. 8, ";")') 

RIN01740 

+RMUI ,NI 



RIN01750 

WRITE(6,'(" ")’) 



RIN01760 

WRITE(6,'(" POI 

=' ' ,E15.8, ' ' ; PFI 

= " ,E15. 8, ";")') 

RIN01770 

+P0I ,PFI 



RIN01780 

WRITE(6,'(" EMOD 

= ",E15.8,"; POIS 

= ",E15. 8, ";")’) 

RIN01790 

+EMOD,POIS 



RIN01800 

WRITE(6,'(1X,72(" 

*"))') 


RIN01810 

U0=2.*NI*SI 



RIN01820 

CONl=6 .*RMUI*UO*ELI /Cl** 2 


RIN01830 

ELOS=ELI/SI 



RIN01840 

AL=(TI*(RI+TI*. 5 ) ) 

**2/12./ EL 1**4 / ( 1 . - POI S**2 ) 

RIN01850 

BET=C0NI*(RI+TI*.5)**2/CI/TI/EM0D 


RIN01860 

eps=ei/eli 



RIN01870 

EL1-EL1I/ELI 



RIN01880 

Cl=Cll/CI*ELI 



RIN01890 

PO=POI/CON1 



RIN01900 

PF=PFI/C0N1 



RIN01910 

END IF 



RIN01920 

IF(ELOS.GT.l.D-5)XLAM 

=.739/ELOS**.585 


RIN01930 

STARV=ELOS*XLAM 



RIN01940 

WRITE(6,'(" IDIM =' 

' ,12, ' ' ; ' ' )' )IDIM 


RIN01950 

WRITE(6,'(" AL ■' 

' ,E15 .8, ' ' ; BET =", 

E15.8,";")' )AL,BET 

RIN01960 

WRITE(6,'(" EPS =' 

',E15.8,"; P0 = ", 

E15.8,";")' )EPS ,P0 

RIN01970 

WRITE(6,'(" ELI =' 

' ,E15.8, ' ' ; PF =", 

E15.8, ' ' ; ' ' )' )EL1,PF 

RIN01980 

WRITE(6,'(" STARV =' 

' ,E15.8, ' ' ; cl =", 

E15.8, ' ' ; ' ' )' )STARV, 

31 RIN01990 

WRITE(6,'(" ")') 



RIN02000 

WRITE(6, ' ( " NDX =' 

1 I5 II II I5 It It I5 
9 9 f f 9 9 9 

,";")’ )NDX 

RIN02010 

WRITE(6,’(" NDP =' 

! 15 » » 1 1 xc 1 ’ f f TC 

9 LJ 9 9 9 9 9 9 

,";")' )NDP 

RIN02020 

WRITE(6,'(" ")') 



RIN02030 

WRITE(6,'(" NMAX =' 

',15,"; NSIG =",I5, 

" ; ' ' )' )NMAX,NSIG 

RIN02040 

WRITE(6 , ' ( " IPR =' 

',15,"; ICAV =" ,15 , 

" ; ' ' )' )IPR,ICAV 

RIN02050 

WRITE(6, ' ( " XF =' 

' , 2 El 5.5) ' )XF(1) ,XF(2) 


RIN02060 

WRITE(6,'(" XB =' 

’ , 2E15 . 5 ) ' )XB( 1 ) ,XB( 2 ) 


RIN02070 

WRITE(6, * ( IX, 72( ' '*' ' 

),/)’) 


RIN02080 

WRITE(6,'(" OUTPUTS' 

’)') 


RIN02090 

WRITE(6, ' (1X,72( ' '*' ' 

))') 


RIN02100 

C 



RIN02110 

IFLAG=0 



RIN02120 

IF(0. .LT.EPS.AND.EPS. 

LT. 1 . )THEN 


RIN02130 

IFLAG=1 



RIN02140 

DX(I)=ELl/FLOAT(NDX(l)) 


RIN02150 

DP( 1 )=ELl/FLOAT(NDP( 1 ) ) 


RIN02160 

DX( 2 ) =( 1 . -EPS ) / FLOAT( NDX( 2 ) ) 


RIN02170 

DP( 2 ) =( 1 . -EPS )/FLOAT(NDP(2 ) ) 


RIN02180 

DX(3)=EPS/FLOAT(NDX(3) ) 


RIN02190 

DP(3)=EPS/FLOAT(NDP(3) ) 


RIN02200 
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non non non 


FILE: RING FORTRAN L4 MTI 


MON 01/20/86 14:39:26 


PAGE 5 OF 38 


C 


c 


ELSE IF(EPS.EQ.1.)THEN 

RIN02210 

IFLAG=2 

RIN02220 

DX(1)=EL1/FL0AT(NDX(1)) 

RIN02230 

DP(1)=EL1/FL0AT(NDP( 1) ) 

RIN02240 

DX(2)=l./FLOAT(NDX(2)) 

RIN02250 

DP(2)=l./FLOAT(NDP(2)) 

RIN02260 

DX(3)=0. 

RIN02270 

DP(3)=0. 

RIN02280 

ELSE IF( 1 . . LT . EPS . AND . EPS . LT . EL1+ 1 )THEN 

RIN02290 

IFLAG=3 

RIN02300 

DX(1)=(EL1+1.-EPS)/FL0AT(NDX(1)) 

RIN02310 

DP(l)=(ELl+l.-EPS)/FLOAT(NDP(l)) 

RIN02320 

DX(2 )=( EPS-1 . )/FLOAT( NDX( 2 ) ) 

RIN02330 

DP(2)=(EPS-l.)/FLOAT(NDP(2)) 

RIN02340 

DX(3)=l./FLOAT(NDX(3)) 

RIN02350 

DP(3)=1 . /FLOAT(NDP(3) ) 

RIN02360 

ELSE IF(EPS.EQ.1.+EL1)THEN 

RIN02370 

IFLAG=4 

RIN02380 

DX(l)=ELl/FLOAT(NDX(l)) 

RIN02390 

DP(l)=ELl/FLOAT(NDP(l)) 

RIN02400 

DX(2)=l./FLOAT(NDX(2)) 

RIN02410 

DP(2)=l./FLOAT(NDP(2)) 

RIN02420 

DX(3)=0. 

RIN02430 

DP(3)=0. 

RIN02440 

END IF 

RIN02450 

IF(IFLAG.EQ.O)STOP 

RIN02460 

RIN02470 

CALL CHECK( AL , BET , EPS , ELI , NDX , NDP , IRON , RECALC ) 

RIN02480 

IF ( RECALC )THEN 

RIN02490 

IF( IFLAG . LE . 2 )THEN 

RIN02500 

CALL CALCD 

RIN02510 

ELSE IF(IFLAG.EQ.3)THEN 

RIN02520 

CALL CALCD3 

RIN02530 

ELSE IF(IFLAG.EQ.4)THEN 

RIN02540 

CALL CALCD4 

RIN02550 

END IF 

RIN02560 

ELSE 

RIN02570 

WRITE(6,'(20X," * NOTE: C,D NOT RECALCULATED ! *”)') 

RIN02580 

END IF 

RIN02590 

RIN02600 

ELASTIC ANALYSIS 

RIN02610 

RIN02620 

WRITE(6,'(" ” ,72( ) 

RIN02630 

WRITE(6,'(26X," * ELASTIC SOLUTION *",/)') 

RIN02640 

CALL ELAS 

RIN02650 

RIN02660 

STEADY STATE ANALYSIS 

RIN02670 

RIN02680 

IELAS=0 

RIN02690 

RIN02700 

FORWARD FLOW SOLUTION 

RIN02710 

RIN02720 

IBACK=0 

RIN02730 

WRITE(6,'(" ”,72("-”))’) 

RIN02740 

WRITE(6,'(23X," * PUMPING FLOW SOLUTION *",/)') 

RIN02750 
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C 

3000 


IF(XF( 1 )**2+XF( 2 )**2 .LT.’l .D-5 )THEN 
XX( 1 )=1 . -HELAS+ .01 
XX(2)=C1-H1ELAS+.01 
ELSE 

XX(1)=1.-XF(1) 

XX(2)=C1-XF(2) 

END IF 

IF(IFLAG.LE.2)THEN 

CALL ZSCNT(EVAL,NSIG,2,NMAX,PAR,XX,FN0RM,WK,IER) 

ELSE IF ( IFLAG . EQ . 3 )THEN 

CALL ZSCNT ( EVAL3 , NS IG , 2 , NMAX , PAR , XX , FNORM , WK , I ER ) 
ELSE IF(IFLAG.EQ.4)THEN 

CALL ZSCNT ( EVAL4 ,NSIG,2, NMAX , P AR , XX , FNORM , WK , I ER ) 

END IF 

IF(IER.NE.O)CALL ERRMSG 
CALL PRTOUT(XX) 

XF(1)=1.-XX(1) 

XF(2)=C1-XX(2) 

FFLOW=RK 

BACK FLOW SOLUTION 
IBACK=1 

WRITE(6 , * ( * * ' ' ,72( * * — ' *))*) 

WRITE(6,'(25X," * BACK FLOW SOLUTION *",/)') 
IF(XB(1)**2+XB(2)**2.LT.1.D-5)THEN 
XX( 1 )=1 . -HELAS 
XX(2)=C1-H1ELAS 
ELSE 

XX(1)=1.-XB(1) 

XX(2)=C1-XB(2) 

END IF 
BFLOW=0. 

I F ( FORCE . LT . 1 . D-8 ) THEN 
IF(IFLAG.LE.2)THEN . 

CALL ZSCNT(EVAL,NSIG, 2, NMAX, PAR.XX, FNORM, WK,IER) 
ELSE I F ( I FLAG . EQ . 3 )THEN 

CALL ZSCNT( EVAL3 , NSIG , 2 , NMAX , PAR , XX , FNORM , WK , IER ) 
ELSE I F ( I FLAG . EQ . 4 )THEN 

CALL Z SCNT ( EVAL4 , NS I G , 2 , NMAX , P AR , XX , FNORM , WK , I ER ) 
END IF 

IF(IER.NE.O)CALL ERRMSG 
BFLOW=RK 
END IF 

CALL PRTOUT(XX) 

XB(1)=1.-XX(1) 

XB(2)=C1-XX(2) 

WRITE(6, , (” ’ ’ ,72( 

CONTINUE 

NFLOW=FFLOW-BFLOW 
WRITE(6,'(32X,"* FLOWS *",/)') 

WRITE( 6, '(5X,' ’PUMPING' ' ,10X,' 'BACK' ' ,11X,"XCAV' ' ,9X 
+ , "NET FLOW", 3X,” STARVED NET FLOW")') 

SFLOW=NFLOW 
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HSTV=XB(1) 

DHSTV=XB(2) 

IF( IBACK . EQ . 1 . AND . FORCE . GT . 1 . D-8 )THEN 
XCAV=1 . . 

HSTV=HELAS 
DHSTV=H1ELAS 
END IF 

IF(XCAV . GT . 1 . D-6 . AND .XCAV . LT . 1 . 00000 1 ) SFLOW=FFLOW*( 1 . + 

+2 .*ELOS*XLAM*XCAV**2*DHSTV/( 2 .*HSTV-DHSTV*( 2 .-XCAV) ) )-BFLOW 
WRITE(6, ' (5E15.6) ' )FFLOW,BFLOW 
♦ ,XCAV , NFLOW , SFLOW 

WRITE(6,'(1X,72(’ '*"),/)' ) 

GOTO 2001 
2000 CALL EXIT 
END 


RIN03310 

RIN03320 

RIN03330 

RIN03340 

RIN03350 

RIN03360 

RIN03370 

RIN03380 
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C> 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c> 


SUBROUTINE AIN(A,B) 


FUNCTION - CALCULATE INVERSE OF 2X2 MATRIX 

RESTRICTIONS 

REMARKS 

EXTERNAL REFERENCES - NONE 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

A INPUT MATRIX 

B OUTPUT MATRIX (INVERSE OF A) 


IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION A(2,2),B(2,2) 

D=A(1,1)*A(2,2)-A(2,1)*A(1,2) 

B(1,1)=A(2,2)/D 

B(2,2)=A(1,1)/D 

B( 1 ,2)=-A( 1 ,2)/D 

B(2,l)=-A(2,l)/D 

RETURN 

END 
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SUBROUTINE AMU(A,B,C) 


FUNCTION 


- PERFORM MATRIX MULTIPICATION OF 


RESTRICTIONS 

REMARKS 

EXTERNAL REFERENCES - NONE 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

A INPUT MATRIX 

B INPUT MATRIX 

C OUTPUT MATRIX (C=AXB) 


IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION A(2,2),B(2,2),C(2,2) 

C( 1 , 1 )=A( l,l)*B(l f l )+A( 1,2)*B(2,1) 

C(2 , 1)=A(2 , I )*B( 1 , 1 )+A(2 ,2)*B(2 , 1) 

C( 1 ,2)=A( 1, 1 )*B( 1 ,2 )+A( I , 2)*B(2 ,2) 

C(2,2)=A(2,1)*B(1,2)+A(2,2)*B(2,2) 

RETURN 

END 
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C> 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C> 


C 


100 


SUBROUTINE CALCD 


FUNCTION - CALCULATE INFLUENCE COEFFICIENTS C,D 

FOR H AND H' AT PSI=1 

RESTRICTIONS - FOR MULTIPLE RUNS, ONLY COMPUTED WHEN 

CERTAIN PARAMETERS CHANGE. ALWAYS 
COMPUTED AT LEAST ONCE. 


REMARKS - NOTE VARIABLES PASSED IN COMMON 

EXTERNAL REFERENCES - AIN,AMU,RUK 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/BCD /XLOC( 100 ) ,C( 100 ) ,D( 100 ) ,NDX(3) ,NDP(3 ) , DX(3 ) ,DP(3 ) 

COMMON /BELAS / FORCE , HELAS , HI ELAS , H2ELAS , H3 ELAS ,W1(4),W2(4) 
COMMON/BPAR / AL , BET , PO , EPS , PF , S , U , DT 
COMMON /BFLAG /IFLAG,IPL0AD 
DIMENSION YO(4) ,YT(4) ,DJ(4) ,CKJ(4,4) 

DIMENSION XLEN(3) 

DIMENSION Al( 2, 2,100 ),A2( 2, 2, 100) 

DIMENSION B1(2,2,100),B2(2,2,10Q) 

DIMENSION A1I(2,2),AT(2,2) 

EXTERNAL DFN1 

XLEN( 1 )=DX( 1 )*NDX( 1 ) 

XLEN(2)=DX(2)*NDX(2) 

XLEN(3)=DX(3)*NDX(3) 

ELl=XLEN( 1 ) 

Y0(l)=0. 

YO(2 )=0. 

YO(3)=l. 

Y0(4)=0. 

CALL RUK(DX( 1 ) ,XLEN( 1 ) ,-ELl ,XNN, YO, Y0,4,DFN1 ,YT,DJ,CKJ,4) 

CALL RUK( DX(2)/2.,DP(2)/2., XNN , XNN , YO , YO , 4 , DFN1 , YT , DJ , CK J , 4 ) 

XLOC(l)=XNN 

Al(l,l,l)=YO(l) 

Al(2,l,l)=YO(2) 

A2(l,l,l)=YO(3) 

A2(2 , 1 , 1 )=YO(4 ) 

DO 100 1=2 ,NDP(2) 

CALL RUK (DX(2),DP(2), XNN , XNN , YO , YO , 4 , DFN 1 , YT , D J , CK J , 4 ) 

XLOC(I )=XNN 
Al( 1 ,1,1 )=YO( 1 ) 

Al(2,l,I)=YO(2) 

A2( 1 , 1 , I )=YO(3) 

A2(2,1,I)=Y0(4) 

CALL RUK(DX(2)/2.,DP(2)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
IF(IFLAG.EQ.2)GOTO 111 

CALL RUK(DX(3)/2.,DP(3)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
XLOC(NDP(2)+l)=XNN 
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200 
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112 
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Al( 1 , 1 ,NDP(2 )+l )=YO( 1 ) 

RIN04570 

Al(2,l,NDP(2)+l)=YO(2) 

RIN04580 

A2( 1 , 1 ,NDP(2 )+l )=YO(3 ) 

RIN04590 

A2(2,l,NDP(2)+l)=YO(4) 

RIN04600 

DO 101 I=2,NDP(3) 

RIN04610 

CALL RUK(DX(3),DP(3),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) . 

RIN04620 

XLOC ( NDP ( 2 ) + I ) =XNN 

RIN04630 

Al( 1 , 1 ,NDP(2 )+I )=YO( 1 ) 

RIN04640 

Al(2,l,NDP(2)+I)=YO(2) 

RIN04650 

A2(l,l,NDP(2)+I)=YO<3) 

RIN04660 

A2(2 , 1 ,NDP(2 )+I )=YO(4) 

RIN04670 

CALL RUK(DX(3)/2.,DP(3)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

RIN04680 

Wl(l)=YO(l) 

RIN04690 

Wl(2)=YO(2) 

RIN04700 

Wl(3)=YO(3) 

RIN04710 

W1(4)=Y0(4) 

RIN04720 


RIN04730 

YO( 1 )=0. 

RIN04740 

Y0(2)=0. 

RIN04750 

Y0(3)=0. 

RIN04760 

YO(4)=l . 

RIN04770 

CALL RUK(DX(1),XLEN(1),-EL1,XNN,Y0,Y0,4,DFN1,YT,DJ,CKJ,4) 

RIN04780 

CALL RUK(DX(2)/2.,DP(2)/2.,XNN,XNN,YO,YO,4,DFN1,YT,DJ,CKJ,4) 

RIN04790 

Al( 1 , 2 , 1 )=YO( 1 ) 

RIN04800 

Al(2,2 : ,l)=YO(2) 

RIN04810 

A2(l,2,l)=YO(3) 

RIN04820 

A2 (2,2,1) =YO( 4 ) 

RIN04830 

DO 200 I=2,NDP(2) 

RIN04840 

CALL RUK(DX(2),DP(2),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

RIN04850 

A1(1,2,I)=Y0(1) 

RIN04860 

A1(2,2,I)=Y0(2) 

RIN04870 

A2(l,2,l)=YO(3) 

RIN04880 

A2(2,2,I)=Y0(4) 

RIN04890 

CALL RUK(DX(2)/2.,DP(2)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

RIN04900 

I F ( I FLAG . EQ . 2 ) GOTO 112 

RIN04910 

CALL RUK(DX(3)/2.,DP(3)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

RIN04920 

A1 ( 1 ,2 ,NDP(2)+1 )=YO( 1 ) 

RIN04930 

Al(2,2,NDP(2)+l)=YO(2) 

RIN04940 

A2( 1 ,2 ,NDP(2 ) + l )=YO( 3 ) 

RIN04950 

A2(2,2,NDP(2)+l)=YO(4) 

RIN04960 

DO 201 I=2,NDP(3) 

RIN04970 

CALL RUK(DX(3),DP(3),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

RIN04980 

Al( 1 ,2,NDP(2)+I )=YO( 1 ) 

RIN04990 

Al(2,2,NDP(2)+I)=YO(2) 

RIN05000 

A2(l,2,NDP(2)+I)=YO(3> 

RIN05010 

A2(2,2,NDP(2)+I)=YO(4) 

RIN05020 

CALL RUK(DX(3)/2. ,DP(3)/2. ,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

RIN05030 

W2( 1 )=YO( 1 ) 

RIN05040 

W2(2)=YO(2) 

RIN05050 

W2(3)=YO(3) 

RIN05060 

W2(4)=YO(4) i 

RIN05070 


RIN05080 

YO( 1 )=1 . 

RIN05090 

YO(2)=0. 

RIN05100 

Y0(3)=0. 

RIN05110 
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Y0(4)=0. 

XNN=1. 

I F ( I FLAG . EQ . 2 ) GOTO 113 

CALL RUK( -DX( 3 ) / 2 . , -DP( 3 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN 1 , YT , DJ , CKJ , 4 ) 
81(l,l,NDP(2)+NDP(3))=YO(l) 

Bl(2,l,NDP(2)+NDP(3))=YO(2) 

B2Q,l,NDP(2)+NDP(3))=YO(3) 

B2(2 , 1 ,NDP(2 )+NDP(3 ) )=YO(4) 

DO 300 II=2,NDP{3) 

I=NDP(2)+NDP(3)-II+1 

CALL RUK(-DX(3),-DP(3),XNN,XNN,Y0,Y0,4,DFN1,YT,DJ,CKJ,4) 
B1(1,1,I)=Y0(1) 

Bl(2,l,I)=YO(2) 

B2(l,l,X)-YO(3) 

300 B2(2,l,I)=YO(4) 

CALL RUK ( -DX( 3 ) / 2 . , -DP ( 3 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN1 , YT , DJ , CKJ , 4 ) 

113 CALL RUK(-DX(2)/2. ,-DP(2)/2 . ,XNN ,XNN, YO, Y0,4,DFN1 , YT,DJ,CKJ ,4 ) 
Bl( 1 ,1 ,NDP(2) )=YO( 1 ) 

Bl(2,l,NDP(2))=YO(2) 

B2( 1 ,1 ,NDP(2) )=YO(3) 

B2(2 , 1 ,NDP(2) )=YO(4) 

DO 301 II=2,NDP(2) 

I=NDP(2)-II+1 

CALL RUK(-DX(2),-DP(2),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

Bl(l , 1 ,1 )=YO( 1 ) 

Bl(2,l,I)=YO(2) 

B2(l,l,D=YO(3) 

301 B2C2,l,I)=YO(4) 

C - 

YO(1)=0. 

YO(2)=l. 

Y0(3)=0. 

Y0(4)=0. 

XNN=1 . 

I F ( I FLAG . EQ . 2 ) GOTO 114 

CALL RUK(-DX(3)/2.,-DP(3)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,£) 
Bl(l,2,NDP(2)+NDP(3))=YO(l) 

Bl(2,2,NDP(2)+NDP(3))=YO(2) 

B2(l,2,NDP(2)+NDP(3))=YO(3) 

B2(2,2,NDP(2)+NDP(3))=YO(4) 

DO 400 11=2 ,NDP(3) 

I=NDP(2)+NDP( 3)-II+l 

CALL RUK(-DX(3),-DP(3),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

Bl( 1 ,2 ,1 )=YO( 1 ) 

Bl(2,2,I)=YO<2) 

B2( 1 ,2,1 )=YO( 3 ) 

400 B2(2,2,I)=YO(4) 

CALL RUK(-DX(3)/2.,-DP(3)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

114 CALL RUK( -DX( 2 ) / 2 . , -DP( 2 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN 1 , YT , DJ , CKJ , 4 ) 
Bl(l,2,NDP(2))=YO(l) 

Bl(2,2,NDP(2))=YO(2) 

B2(l,2,NDP(2))=YO(3) 

B2(2,2,NDP(2) )=YO(4) 

DO 401 II=2,NDP(2) 

I=NDP(2 )-II+l 
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CALL RUK(-DX(2),-DP(2),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl( 1 ,2,1 )=YO( 1 ) 

B1(2,2,I)=Y0(2) 

B2(1,2,I)=Y0(3) 

401 B2(2,2,I)=Y0(4) 

C 

DO 500 1=1 ,NDP(2 ) 

CALL AIN(A1( 1 ,1,1), All) 

CALL AMU(A2( 1 , 1 , I ) , All ,AT) 

CALL AMU(AT,B1(1,1,I),A1(1,1,I)) 

DO 21 11=1,2 
DO 21 JJ=1 , 2 

21 AT(II , JJ)=B2( II , JJ, I )-Al( II , JJ, I ) 

CALL AIN(AT,B2( 1,1,1)) 

C(I )=B2( 1 ,2,1 ) /AL*DP(2) 

D(I )=B2(2 ,2 , I )/AL*DP(2 ) 

500 CONTINUE 

I F ( I FLAG . EQ . 2 ) RETURN 
INDEX1=NDP(2)+1 
INDEX2=NDP( 2 )+NDP( 3 ) 

DO 501 I=INDEX1 , INDEX2 
CALL AIN(A1( 1 ,1,1), All) 

CALL AMU (A2 (1,1, I), All, AT) 

CALL AMU(AT,B1(1,1,I),A1(1,1,I)) 

DO 22 11=1,2 
DO 22 JJ=1 ,2 

22 AT(II,JJ)=B2(II,JJ,I)-A1(II,JJ,I) 

CALL AIN(AT,B2( 1,1,1)) 

C ( I ) =B2 (1,2,1)/ AL*DP( 3 ) 

D( I )=B2 (2,2,1)/ AL*DP( 3 ) 

501 CONTINUE 
C 

RETURN 

END 
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C> 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C> 


C 


100 


SUBROUTINE CALCD3 


FUNCTION - CALCULATE INFLUENCE COEFFICIENTS C,D 

FOR H AND H' AT PSI=1 

RESTRICTIONS - FOR MULTIPLE RUNS, ONLY COMPUTED WHEN 

CERTAIN PARAMETERS CHANGE. ALWAYS 
COMPUTED AT LEAST ONCE. 


REMARKS - NOTE VARIABLES PASSED IN COMMON 

EXTERNAL REFERENCES - AIN,AMU,RUK 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/BCD /XLOC( 100) ,C( 100) ,D( 100) ,NDX( 3 ) ,NDP(3 ) ,DX( 3 ) ,DP(3 ) 
COMMON/ BELAS / FORCE , HELAS , HI ELAS , H2 ELAS , H3 ELAS , W1 ( 4 ) , W2 ( 4 ) 
COMMON/ BPAR /AL,BET,PO,EPS, PF,S,U,DT 
COMMON /B FLAG /lFLAG,IPLOAD 
DIMENSION YO(4),YT(4),DJ(4),CKJ(4,4) 

DIMENSION XLEN(3) 

DIMENSION Al(2 ,2 , 100 ) ,A2(2 ,2 , 100) 

DIMENSION B1(2,2,100),B2(2,2,100) 

DIMENSION A1I(2,2),AT(2,2) 

EXTERNAL DFN1 

XLEN( 1 )=DX( 1 )*NDX( 1 ) 

XLEN(2)=DX(2)*NDX(2) 

XLEN( 3 )=DX( 3 )*NDX( 3 ) 

ELl=XLEN(l)+XLEN(2) 

YO(1)=0. 

Y0(2)=0. 

YO(3)=l . 

YO(4)=0 . 

CALL RUK(DX( 1 ) ,XLEN( 1 ) ,-ELl ,XNN, YO,YO,4 ,DFN1 , YT,DJ,CKJ ,4) 

CALL RUK( DX( 2)/2.,DP(2)/2., XNN , XNN , YO , YO , 4 , DFN1 , YT , DJ , CKJ , 4 ) 
XLOC( I)=XNN 
Al(l ,1,1 )~YO( 1 ) 

Al(2,l,l)=YO(2) 

A2(l,l,l)=YO(3) 

A2(2 , 1 , 1 )==YQ(4) 

DO 100 I=2,NDP(2) 

CALL RUK(DX( 2 ) ,DP(2 ) ,XNN,XNN,YO,YO,4 , DFN1 ,YT,DJ,CKJ ,4) 

XLOC(l)=XNN 

A1(1,1,I)=Y0(1) 

A1(2,1,I)*Y0(2) 

A2(l,l,I)=YO(3) 

A2(2,l,I)*YO(4) 

CALL RUK (DX(2)/2.,DP(2)/2., XNN , XNN , YO , YO , 4 , DFN 1 , YT , D J , CKJ , 4 ) 
CALL RUK( DX( 3)/2.,DP(3)/2., XNN , XNN , YO , YO , 4 , DFN1 , YT , DJ , CKJ , 4 ) 
XLOC(NDP(2)+l)=XNN 
Al( 1 , 1,NDP(2 )+l )=YO( 1) 
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Al(2,l,NDP(2)+l)=YO(2) 

RIN06560 

A2( 1 , 1 ,NDP(2)+1 )=YO(3) 
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A2(2,l,NDP(2)+l)=YO(4) , 
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DO 101 I=2,NDP(3) 
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CALL RUK(DX(3),DP(3),XNN,XNN,Y0,Y0,4,DFN1,YT,DJ,CKJ,4) 
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RIN06870 

CALL RUK(DX(2)/2.,DP(2)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

RIN06880 
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RIN06950 

Al(l ,2,NDP(2)+I )=YO( 1 ) 

RIN06960 

Al(2,2,NDP(2)+I)=YO(2) 

RIN06970 

A2(l,2,NDP(2)+I)=YO(3) 

RIN06980 

A2(2,2,NDP(2)+I)=YO(4) 

RIN06990 

CALL RUK(DX(3)/2. ,DP(3)/2. ,XNN,XNN,Y0,Y0,4,DFN1,YT,DJ,CKJ,4) 

RIN07000 

W2( 1 )=YO( 1 ) 

RIN07010 

W2(2)=YO(2) 

RIN07020 

W2(3)=YO(3) 

RIN07030 

W2(4)=YO(4) 

RIN07040 


RIN07050 

YO(l)=l. 

RIN07060 

Y0(2)=0. 

RIN07070 

YO(3)=0. 

RIN07080 

YO(4)=0 . 

RIN07090 

XNN=1. 

RIN07100 


142 



PAGE 15 OF 38 


FILE: RING FORTRAN L4 MTI MON 01/20/86 14:39:26 


300 


301 

C 


400 

114 


CALL RUK(-DX(3)/2.,-DP(3)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl(l,l,NDP(2)+NDP(3))=YO(l) 

Bl(2,l,NDP(2)+NDP(3))=YO(2) 

B2(l,l,NDP(2)+NDP(3))=YO(3) 

B2(2,l,NDP(2)+NDP(3))=YO(4) 

DO 300 II=2,NDP(3) 

I=NDP(2)+NDP(3)-II+1 

CALL RUK( -DX( 3 ) , -DP( 3 ) , XNN , XNN , YO , YO , 4 , DFN1 , YT , D J , CK J , 4 ) 
Bl(l,l,I)=YO(l) 

Bl(2,l,I)=YO(2) 

B2(l,l,I)=YO(3) 

B2(2,l,I)=YO(4) 

CALL RUK ( -DX( 3 ) / 2 . , -DP( 3 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN 1 , YT , DJ , CK J , 4 ) 
CALL RUK ( -DX ( 2 ) / 2 . , -DP ( 2 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN 1 , YT , D J , CK J , 4 ) 
Bl(l,l,NDP(2))=YO(l) 

Bl(2,l,NDP(2))=YO(2) 

B2( 1 , 1 ,NDP(2 ) )=YO(3 ) 

B2(2 , 1 ,NDP(2 ) )=Y0(4) 

DO 301 II=2,NDP(2) 

I=NDP(2)-II+1 

CALL RUK(-DX(2),-DP(2),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl(l,l,I)=YO(l) 

Bl(2,l,I)=YO(2) 

B2(l,l,I)=YO(3) 

B2(2 f l,I)=YO(4) 

YO(1)=0. 

YO(2)=l. 

Y0(3)=0. 

YO(4)=0. 

XNN=1. 

CALL RUK ( — DXC 3 ) / 2 . , -DP ( 3 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN 1 , YT , D J , CK J , 4 ) 
Bl(l,2,NDP(2)+NDP(3))=YO( 1 ) 

Bl(2,2,NDP(2)+NDP(3))=YO(2) 

B2(l,2,NDP(2)+NDP(3))=YO(3) 

B2(2,2,NDP(2)+NDP(3))=YO(4) 

DO 400 11=2 ,NDP( 3) 

I=NDP(2)+NDP(3)-II+1 

CALL RUK(-DX(3),-DP(3),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl(l,2,I)=YO(l) 

Bl(2,2,I)=YO(2) v 

B2(l,2,I)=YO(3) ‘ 

B2(2,2,I)=YO(4) =: 

CALL RUK( -DX (3)/2.,-DP(3)/2., XNN , XNN , YO , YO , 4 , DFN 1 , YT , D J , CK J , 4 ) 
CALL RUK(-DX(2)/2.,-DP(2)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl( 1 ,2,NDP(2 ) )=YO( 1 ) 

Bl(2,2,NDP(2))=YO(2) 

B2(l,2,NDP(2))=YO(3) 

B2 ( 2 , 2 , NDP ( 2 ) ) =YO( 4 ) 

DO 401 II=2,NDP(2) 

I=NDP(2)-II+1 

CALL RUK(-DX(2),-DP(2),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl(l,2,I)=YO(l) 

Bl(2,2,I)=YO(2) 

B2(l,2,I)=YO(3) 
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401 B2(2,2,I)=YO(4) 

C 

DO 500 1=1 ,NDP(2) 

CALL AIN(A1(1,1,I),A1I) 

CALL AMU(A2( 1 ,1,1), All, AT) 

CALL AMU(AT,B1 (1, 1,1 ),A1 (1,1,1)) 
DO 21 11=1,2 
DO 21 JJ=1 ,2 

21 AT(II,JJ)=B2(II,JJ,I)-A1(II,JJ,I) 
CALL AIN(AT,B2( 1,1,1)) 
C(I)=B2(1,2,I)/AL*DP(2) 
D(I)=B2(2,2,I)/AL*DP(2) 

500 CONTINUE 
INDEXl=NDP(2)+l 
INDEX2=NDP(2 )+NDP( 3 ) 

DO 501 I=INDEX1 , INDEX2 
CALL AIN(A1( 1 ,1,1), All) 

CALL AMU(A2( 1 ,1,1), All, AT) 

CALL AMU(AT,B1(1,1,I),A1(1,1,I)) 
DO 22 11=1,2 
DO 22 JJ=1 ,2 

22 AT(II,JJ)=B2(II,JJ,I)-A1(II,JJ,I) 
CALL AIN(AT,B2( 1 ,1,1)) 

C( I )=B2 (1,2,1)/ AL*DP( 3 ) 

D( I )=B2 (2,2,1)/ AL*DP( 3 ) 

501 CONTINUE 
C 

RETURN 

END 
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100 


SUBROUTINE CALCD4 


FUNCTION 


CALCULATE INFLUENCE COEFFICIENTS C,D 
FOR H AND H' AT PSI=1 


RESTRICTIONS - FOR MULTIPLE RUNS, ONLY COMPUTED WHEN 

CERTAIN PARAMETERS CHANGE. ALWAYS 
COMPUTED AT LEAST ONCE. 


REMARKS 


NOTE VARIABLES PASSED IN COMMON 


EXTERNAL REFERENCES - AIN,AMU,RUK 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/BCD /XLOC( 100 ) , C( 100 ) , D< 100 ) ,NDX( 3 ) ,NDP( 3 ) ,DX( 3 ) , DP( 3 ) 

COMMON/BELAS /FORCE, HELAS , H1ELAS , H2ELAS ,H3ELAS , W1 (4 ) , W2(4 ) 
COMMON/BPAR /AL,BET,PO,EPS,PF,S,U,DT 
COMMON /BFLAG /I FLAG, IPLOAD 
DIMENSION YO(4) ,YT(4) ,DJ(4) ,CKJ(4,4) 

DIMENSION XLEN(3) 

DIMENSION A1(2,2,100),A2(2,2,100) 

DIMENSION B1(2,2,100),B2(2,2,100) 

DIMENSION A1I(2,2),AT(2,2) 

EXTERNAL DFN1 

XLEN( 1 )=DX( 1 )*NDX( 1 ) 

XLEN(2)=DX(2)*NDX(2) 

XLEN(3)=DX(3)*NDX(3) 

EL1=XLEN( 1 ) 

YO( 1 )=0. 

YO(2)=0. 

YO(3)=l. 

Y0(4)=0. 

CALL RUK(DX(l)/2. ,DP( 1 )/2 . ,-ELl ,XNN, YO, YO,4 ,DFN1 , YT,DJ,CKJ,4) 
XLOC( 1 )=XNN 
A1 ( 1 , 1 , 1 )=YO( 1 ) 

Al(2,l,l)=YO(2) 

A2( 1,1,1 )=YO(3) 

A2(2 , 1 , 1 )=YO(4) 

DO 100 1=2 ,NDP( 1 ) 

CALL RUK(DX( 1 ) ,DP(1 ) ,XNN,XNN, YO, YO,4 ,DFN1 , YT,DJ,CKJ,4) 

XLOC(I)=XNN 

Al( 1 , 1 , 1 )=YO( 1 ) 

Al(2 , 1,1 )=YO(2 ) 

A2(l,l,I)=YO(3) 

A2(2, l,I)=YO(4) 

CALL RUK ( DX ( 1 ) / 2 . , DP ( 1 ) / 2 . , XNN , XNN ,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
CALL RUK (DX(2)/2.,DP(2)/2., XNN , XNN , YO , YO , 4 , DFN 1 , YT , D J , CK J , 4 ) 
XLOC(NDP{ 1 )+l )=XNN 
Al( 1 , 1 ,NDP( 1 )+l )=YO( 1 ) 

Al(2 , 1 ,NDP( 1 )+l )=YO(2 ) 
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200 


201 
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A2(l,l,NDP(l)+l)=YO(3) 

A2(2, i,NDP( 1 ) + l )=Y0(4) 

DO 101 1=2 f NDP(2 ) 

CALL RUK(DX(2).,DP<2),XNN,XNN,Y0,Y0,4,DFN1,YT,DJ,CKJ,4) 

XL0C(NDP(1)+I)=XNN 

Al( 1 , 1 ,NDP( 1 ) + I )=Y0( 1) 

A1 (2 , 1 ,NDP( 1 ) + I )=YO(2 ) 

A2( 1 , 1 ,NDP( 1 ) + I )=YO(3) 

A2(2 , 1 ,NDP( 1 ) + I )=Y0(4) 

CALL RUK(DX(2)/2.,DP(2)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

W1 ( 1 )=YO( 1 ) 

Wl(2)=YO(2) - 

Wl(3)=YO(3) 

Wl(4)=YO(4) 

YO( 1 )=0 . 

Y0(2)=0. 

YO(3)=0. 

YO(4)=l . 

CALL RUK(DX(l)/2.,DP(l)/2.,-ELl,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Al(l,2,l)=YO(l) 

Al(2,2,l)=YO(2) 

A2(l,2,l)=YO(3) 

A2(2,2,l)=YO(4) 

DO 200 1=2 ,NDP( 1 ) 

CALL RUK(DX(l),DP(l),XNN f XNN,Y0,Y0,4,DFNl,YT,DJ,CKJ,4) 
Al(l,2,I)=YO(l) 

Al(2,2,I)=YO(2) 

A2( 1 ,2,1 )=YO(3) 

A2(2,2,I)=YO(4) 

CALL RUK(DX(l)/2.,DP(l)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
CALL RUK (DX(2)/2.,DP(2)/2., XNN , XNN , YO , YO , 4 , DFN 1 , YT , D J , CK J , 4 ) 
Al(l,2,NDP(l)+l)=YO(l) 

Al(2,2,NDP(l)+l)=YO(2) 

A2( 1 ,2 ,NDP( 1 )+l )=YO( 3 ) 

A2(2 ,2 ,NDP( 1 ) + l )=YO(4) 

DO 201 I=2,NDP(2) 

CALL RUK(DX(2),DP(2),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

Al( 1 ,2,NDP( 1 )+I )=YO( 1 ) 

Al(2,2,NDP(l)+I)=YO(2) 

A2(l f 2,NDP(l)+I)=YO(3) 

A2( 2 ,2 ,NDP( 1 )+I )=Y0(4) 

CALL RUK( DX( 2)/2.,DP(2)/2., XNN , XNN , YO , YO , 4, DFN 1 , YT , DJ , CKJ , 4 ) 
W2( 1 )=YO( 1 ) 

W2(2)=YO(2) 

W2(3)=YO(3) 

W2(4)=YO(4) 

YO(l)=l. 

YO(2)=0. 

YO(3)=0. 

Y0(4)=0. 

XNN=1. • 

CALL RUK(-DX(2)/2.,-DP(2)/2.,XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl(l,l,NDP(l)+NDP(2))=YO(l) 
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300 


301 

C 


400 


401 

C 


Bl(2,l,NDP(l)+NDP(2))=YO(2) 

B2( 1 , 1 ,NDP( 1 )+NDP(2) )=YO(3) 

B2 (2,1 ,NDP( 1 )+NDP(2 ) )=YO(4) 

DO 300 II=2,NDP(2) 

I=NDP(1)+NDP(2)-II+1 

CALL RUK(-DX(2),-DP(2),XNN,XNN,Y0,Y0,4,DFN1,YT,DJ,CKJ,4) 
Bl(l,l,I)=YO(l) 

Bl(2,l,I)=YO<2) 

B2(l,l,I)=YO(3) 

B2(2,l,I)=YO(4) 

CALL RUK( -DX( 2 ) / 2 . , -DP( 2 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN1 , YT , DJ , CK J , 4 ) 
CALL RUK ( -DX( 1 ) / 2 . , -DP( 1 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN1 , YT , DJ , CK J , 4 ) 
B 1 ( 1 , 1 ,NDP( 1 ) )=YO( 1 ) 

B1(2,1,NDP( 1 ) )=YO(2) 

B2(l,l,NDP(l))=YO(3) 

B2(2 , 1 ,NDP( 1 ) )=YO(4) 

DO 301 11=2 ,NDP( 1) 

I=NDP( 1 )-II+l 

CALL RUK(-DX(1),-DP(1),XNN,XNN,Y0,Y0,4,DFN1,YT,DJ,CKJ,4) 

Bl(l, l,I)=YO( 1 ) 

Bl(2,l,I)=YO(2) 

B2(l,l,I)=YO(3) 

B2(2,l,I)=YO(4) 

YO(1)=0. 

YO(2)=l. 

Y0(3)=0. 

YO(4)=0. 

XNN=1. 

CALL RUK (-DX(2)/2.,-DP(2)/2., XNN , XNN f YO , YO , 4 , DFN 1 , YT , D J , CK J , 4 ) 
Bl(l,2,NDP(l)+NDP(2))=YO(l) 

Bl(2,2,NDP(l)+NDP(2))=YO(2) 

B2( 1 ,2 ,NDP( 1 )+NDP(2 ) )=YO(3 ) 

B2(2,2 ,NDP( 1 )+NDP(2) )=YO(4) 

DO 400 11=2 ,NDP(2 ) 

I=NDP( 1 )+NDP(2 )-II+l 

CALL RUK(-DX(2),-DP(2),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 
Bl(l,2,I)=YO(l) 

Bl(2,2,I)=YO(2) 

B2(l,2,I)=YO(3) 

B2(2,2,I)=Y0(4) 

CALL RUK ( -DX( 2 ) / 2 . , -DP ( 2 ) / 2 . , XNN , XNN , YO , YO , 4 , DFN 1 , YT , DJ , CK J , 4 ) 
CALL RUK(-DX(l)/2.,-DP(l)/2. ,XNN,XNN, YO, YO,4,DFNl ,YT,DJ,CKJ,4) 
Bl( 1 ,2 ,NDP( 1 ) )=YO( 1 ) 

Bl(2,2,NDP(l))=YO(2) 

B2( 1 ,2 ,NDP( 1 ) )=YO(3) 

B2(2 ,2 ,NDP( 1 ) )=YO(4) 

DO 401 11=2 ,NDP( 1 ) 

I=NDP( 1 )-II+l 

CALL RUK(-DX(l),-DP(l),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

Bl( 1 ,2 ,1 )=YO( 1 ) 

Bl(2,2,I)=YO(2) 

B2(l,2,I)=YO(3) 

B2(2,2,I)=YO(4) 


RIN09050 
RIN09060 
RIN09070 
RIN09080 
RIN09090 
RIN09100 
RIN09110 
RIN09120 
RIN09130 
RIN09140 
RIN09150 
RIN09160 
RIN09170 
RIN09180 
RIN09190 
RIN09200 
RIN09210 
RIN09220 
RIN09230 
RIN09240 
RIN09250 
RIN09260 
RIN09270 
RIN09280 
RIN09290 
RIN09300 
RIN09310 
RIN09320 
RIN09330 
RIN09340 
RIN09350 
RIN09360 
RIN09370 
RIN09380 
RIN09390 
RIN09400 
RIN09410 
RIN09420 
RIN09430 
RIN09440 
RIN09450 
RIN09460 
RIN09470 
RIN09480 
RIN09490 
RIN09500 
RIN095 10 
RIN09520 
RIN09530 
RIN09540 
RIN09550 
RIN09560 
RIN09570 
RIN09580 
RIN09590 


FILE: RING FORTRAN L4 MTI 


MON 01/20/86 14:39:26 


PAGE 20 OF 38 


DO 500 I=1,NDP(1) 

CALL AIN(A1(1,1,I),A1I) 

CALL AMU(A2( 1 ,1,1), All, AT) 

CALL AMU(AT,B1(1,1,I),A1(1,1,I)) 
DO 21 11=1,2 
DO 21 JJ=1 , 2 

21 AT(II , JJ)=B2(II , JJ, I )-Al( II , JJ, I ) 
CALL AIN(AT,B2(1 , 1 , I) ) 

C(I )=B2( 1 ,2 , I ) /AL*DP( 1 ) 

D( I ) =B2 (2,2,1) / AL*DP( 1 ) 

500 CONTINUE 
INDEX1=NDP( 1 )+l 
INDEX2=NDP( 1 )+NDP(2 ) 

DO 501 I=INDEX1 , INDEX2 
CALL AIN(A1(1 , 1 , I ) ,AlI) 

CALL AMU(A2(1 ,1,1), All, AT) 

CALL AMU(AT,B1(1,1,I),A1(1,1,I)) 
DO 22 11=1,2 
DO 22 JJ=1,2 

22 AT( II , JJ)=B2( II , JJ, I )-Al( II , JJ, I ) 
CALL AIN(AT,B2( 1 ,1,1)) 

C( I )=B2( 1 ,2,1) / AL*DP( 2 ) 

D(I )=B2(2 ,2 , I )/AL*DP(2 ) 

501 CONTINUE 
C 

RETURN 

END 


RIN09600 

RIN09610 

RIN09620 

RIN09630 

RIN09640 

RIN09650 

RIN09660 

RIN09670 

RIN09680 

RIN09690 

RIN09700 

RIN09710 

RIN09720 

RIN09730 

RIN09740 

RIN09750 

RIN09760 

RIN09770 

RIN09780 

RIN09790 

RIN09800 

RIN09810 

RIN09820 

RIN09830 

RIN09840 

RIN09850 

RIN09860 
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SUBROUTINE CHECK(AL, BET, EPS , ELI ,NDX,NDP, IRUN, RECALC) RIN09870 

C> RIN09880 

C RIN09890 

C RIN09900 

C FUNCTION - CHECK INPUT VARIABLES FOR RECALCULATION RIN09910 

C INFLUENCE COEFFICIENTS C,D RIN09920 

C RIN09930 

C RESTRICTIONS - ONLY USED FOR MULTIPLE RUNS WITH IRUN.GE.2 RIN09940 

C RIN09950 

C REMARKS - RIN09960 

C RIN09970 

C EXTERNAL REFERENCES - RIN09980 

C RIN09990 

C ARGUMENT DEFINITION: RIN10000 

C NAME DESCRIPTION RIN10010 

C AL,BET,... INPUT VARIABLES, ARE THEY CHANGED ? RIN10020 

C RECALC LOGICAL VARIABLE, IF .TRUE., RECALCULATE C,D RIN10030 

C RIN10040 

C RIN10050 

c > RIN10060 

IMPLICIT REAL*8 (A-H,0-Z) RIN10070 

COMMON/ BOLD /ALOLD, BETOLD , EPSOLD, EL10LD,NDX0LD( 3 ) , NDPOLD( 3 ) RIN 10080 

DIMENSION NDX(3) ,NDP(3) RIN10090 

LOGICAL RECALC, TEMP RIN10100 

RECALC=.TRUE. RIN10110 

IF(IRUN.GT. 1 )THEN RIN10120 

TEMP= . TRUE . RIN10130 

TEMP=TEMP . AND . ( AL .EQ. ALOLD) RIN10140 

TEMP=TEMP. AND. ( EPS. EQ. EPSOLD). AND. ( ELI. EQ.EL10LD) RIN10150 

DO 1 1=1,3 RIN10160 

TEMP=TEMP . AND . ( NDX ( I ) . EQ . NDXOLD (I)) RIN10170 

TEMP=TEMP. AND. (NDP( I ) . EQ . NDPOLD( I ) ) RIN10 180 

1 CONTINUE RIN10190 

IF(TEMP)RECALC=. FALSE. RIN10200 

END IF RIN10210 

ALOLD =AL RIN10220 

EPSOLD=EPS RIN10230 

EL10LD=EL1 RIN10240 

DO 2 1=1,3 RIN10250 

NDXOLD ( I ) =NDX ( I ) ' RIN10260 

NDPOLD (I)=NDP(l) RIN10270 

2 CONTINUE RIN10280 

RETURN RIN10290 

END RIN10300 
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SUBROUTINE 

I CONST(HO,H1,PF) 

RIN10310 




— K1NIU.5ZU 

c 



RIN10330 

c 



RIN10340 

c 

FUNCTION 

- CALCULATE CONSTANTS RK,RC,R1,R2,XCAV,ICAV 

RIN10350 

c 


FOR SLIDER BEARIN PRESSURES 

RIN10360 

c 



RIN10370 

c 

RESTRICTIONS - 

RIN10380 

c 



RIN10390 

c 

REMARKS 

- CALLED BY SUBROUTINE EVAL 

RIN10400 

c 


NOTE OTHER VARIABLES PASSED IN COMMON 

RIN10410 

c 



RIN10420 

c 

EXTERNAL 

REFERENCES - DABS,DSQRT 

RIN10430 

c 



RIN10440 

c 

ARGUMENT 

DEFINITION: 

RIN10450 

c 

NAME 

DESCRIPTION 

RIN10460 

c 

HO 

FILM THICKNESS AT PSI=0 

RIN10470 

c 

HI 

FILM THICKNESS AT PSI=1 

RIN10480 

c 

PF 

RESERVOIR PRESSURE 

RIN10490 

c 



RIN10500 

c 

c> — 



RIN10510 

-RIN10520 


RIN10530 

RIN10540 

RIN10550 

RIN10$60 

RIN10570 

RIN10580 

RIN10590 

RIN10600 

RIN10610 

RIN10620 

RIN10630 

RIN10640 

RIN10650 

RIN10660 

RIN10670 

RIN10680 

RIN10690 

RIN10700 

RIN10710 

RIN10720 

RINI0730 

RIN10740 

RIN10750 

RIN10760 

RINI0770 

RIN10780 

RIN10790 


IMPLICIT REAL*8 (A-H,0-Z) 

* COMMON/ BCOEFF / RK , RC , R 1 , R2 , XCAV , I CAV 
COMMON/BINT / I ELAS , I BACK , I PR 
R0=H1+H0 
R1-H0 
R2=H1-H0 

IF(IBACK.EQ.0)THEN 

RK=2 . *H0*H1/R0*( 1 . -H0*H1*PF ) 

RC=1 . /R0*( 1 . +H1**2*PF ) 

XCAV=0 . 

ELSE I F ( I BACK .EQ.I)THEN 
ALF=Hl*(-R2)*PF 

RK=H1*( 1 .+ALF+DSQRT(DABS(ALF**2+2.*ALF+1 .E-7 ) ) ) 
IF(DABS(R2).GT.1.D-10)THEN 
XCAV=(RK-H1)/R2+1. 

ELSE 

XCAV=1 .0 
END IF 

if((rk.le.ho.and.icav.ne'.o).or.icav.eq.i)then 

RC=l./2./RK 

ELSE 

, RK=2.*H0*H1/R0*(1.+H0*H1*PF) 
RC=1./R0*(1.-HI**2*PF) 

END IF 
END IF 
RETURN 
END 
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C> 

C 

C 

C 

C 

C 

C> 


C 


C 

10 


20 


SUBROUTINE CONST2(HO ,H1,PF) 


REMARKS - SEE CONST 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/BCD /XLOC( 100) ,C(100) ,D( 100 ) ,NDX(3) ,NDP(3 ) ,DX(3 ) ,DP(3) 

COMMON/BCOEFF/RK,RC,Rl,R2,XCAV,ICAV 

COMMON/ BINT /IELAS,IBACK,IPR 

R0=H1+H0 

Rl-HO 

R2=H1-H0 

CTEMP= ( - 1 ) ** ( I BACK+ 1 ) 

RK==2.*HO*H1/RO*( 1 . +CTEMP*H0*H1*PF ) 

RC-1 . /R0*( 1 .-CTEMP*H1**2*PF ) 

XCAV=0 . 

IF(IBACK.EQ. 1)THEN 

IF(IPR.EQ. 1 ) WRITE (4,*) ' ' 

PMIN=100. 

DO 10 1=1 ,NDP(2)+NDP(3) 

IF(IPR.EQ. 1 )WRITE(4,*)P(XL0C(I ) ) 

PMIN=DMIN1 ( PMIN , P(XLOC (I ) ) ) 

CONTINUE 

IF ( PMIN. GE.O. DO) RETURN 

H11=H1-H0 

XCAV=0 . 5 

HCAV=H1+H11*(XCAV-1 . ) 

ICOUNT=0 

ICOUNT=ICOUNT+l 

HCVOLD=HCAV 

FHCAV=PF+1 . /Hll*(-1 . /Hl+HCAV/2 . /Hl**2+1 . / 2. /HCAV) 

F1HCAV=1 . /Hll*( 1 . /2 . /Hl**2-1 . /2 . /HCAV**2 ) 

HCAV=HCAV-FHCAV/ F 1HCAV 
XCAV=(HCAV-H1)/H11+1. 

IF(IPR.EQ. 1 )WRITE(4,*)XCAV 

I F (I COUNT . LT . 1 0 . AND . DABS ( HCAV-HCVOLD ) . GT . 1 . D- 7 )GOTO 2 0 
RK=HCAV 
RC=l./2./HCAV 
END IF 
RETURN 
END 


RIN10800 
-RIN10810 
RIN10820 
RIN10830 . 
RIN10840 
RIN10850 
RIN10860 ,• 
-RIN10870 
RIN10880 
RIN10890 
RIN10900 
RIN10910 . 
RIN10920 ; 
RIN10930 
RIN10940 
RIN10950 
RIN10960 
RIN10970 
RIN10980 
RIN10990 
RIN11000 
RIN11010 
RIN11020 
RIN11030 
RIN11040 
RIN11050 
RIN11060 
RIN11070 
RIN11080 
RIN11090 
RIN11100 
RIN11110 
RIN11120 
RIN11130 
RIN11140 
RIN11150 
RIN11160 
RIN11170 
RIN11180 
RIN11190 
RIN11200 
RIN11210 
RIN11220 
RIN11230 
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C> 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c> 


SUBROUTINE DFNl(X,Y,D) 


FUNCTION - DERIVATIVE ROUTINE USED BY RUK 

(NO HYDRODYNAMICS) 

RESTRICTIONS 

REMARKS - NOTE VARIABLES PASSED IN COMMON 

EXTERNAL REFERENCES - 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

X VALUE OF X FOR DERIVATIVES 

Y VALUES OF INDEPENDENT VARIABLES 

D VALUES OF DERIVATIVES 


IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION Y(4),D(4) 

COMMON/BPAR /AL,BET,PO,EPS,PF,S,U,DT 
D( 1 )=Y(2) 

D(2)=Y(3) 

D(3)=Y(4) 

D(4)=-Y(1)/AL 

RETURN 

END 


RIN11240 

■RIN11250 

RIN11260 

RIN11270 

RIN11280 

RIN11290 

RIN11300 

RIN11310 

RIN11320 

RIN11330 

RIN11340 

RIN11350 

RIN11360 

RIN11370 

RIN11380 

RIN11390 

RIN11400 

RIN11410 

RIN11420 

RIN11430 

RIN11440 

RIN11450 

RIN11460 

RIN11470 

RIN11480 

RIN11490 

RIN11500 

RIN11510 

RIN11520 

RIN11530 
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SUBROUTINE DFN2(X,Y,D) 

c > 

C 

C 

C FUNCTION - DERIVATIVE ROUTINE USED BY RUK 

C (HYDRODYNAMICS INCLUDED) 

C 

C RESTRICTIONS 

C 

C REMARKS - NOTE VARIABLES PASSED IN COMMON 

C 

C EXTERNAL REFERENCES - P 

C 

C ARGUMENT DEFINITION: 

C NAME DESCRIPTION 

C X VALUE OF X FOR DERIVATIVES 

C Y VALUES OF INDEPENDENT VARIABLES 

C D VALUES OF DERIVATIVES 

C 
C 

C> 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION Y(4),D(4) 

COMMON/BPAR / AL , BET , PO , EPS , PF , S , U , DT 
COMMON /BFLAG /IFLAG, IPLOAD 
D( 1 )=Y(2) 

D(2)=Y(3) 

D(3)=Y(4) 

H=0 

IF(IPLOAD.EQ. 1 )H=PO 
IF(X.GT.O.)H=H-P(X) 

H=H*BET 

D(4)=(H-Y(1))/AL 

RETURN 

END 


RIN11540 

■RIN11550 

RIN11560 

RIN11570 

RIN11580 

RIN11590 

RIN11600 

RIN1I610 

RIN11620 

RIN11630 

RIN11640 

RIN11650 

RIN11660 

RIN11670 

RIN11680 

RIN11690 

RIN11700 

RIN117I0 

RIN11720 

RIN11730 

•RIN11740 

RIN11750 

RIN11760 

RIN11770 

RIN11780 

RIN11790 

RIN11800 

RIN11810 

RIN11820 

RIN11830 

RIN11840 

RIN11850 

RIN11860 

RIN11870 

RIN11880 
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C> 

C 

C 

C 

C 

C , 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C>- 


SUBROUTINE DFN3(X,Y,D) 


FUNCTION - DERIVATIVE ROUTINE USED BY RUK 

(HYDRODYNAMICS INCLUDED, NO PRELOAD PO) 

RESTRICTIONS - 

REMARKS - NOTE VARIABLES PASSED IN COMMON 

EXTERNAL REFERENCES - P 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

X VALUE OF X FOR DERIVATIVES 

Y VALUES OF INDEPENDENT VARIABLES 

D VALUES OF DERIVATIVES 


IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION Y(4),D(4) 

COMMON/BPAR /AL,BET,PO,EPS,PF,S,U,DT 
COMMON/BINT /IELAS,IBACK,IPR 
D( 1 )=Y(2 ) 

D(2)=Y(3) 

D(3)=Y(4) 

H=BET*P0 

D(4)=(H-Y( 1 ) )/AL 

RETURN 

END 


RIN11890 

■RIN11900 

RIN11910 

RIN11920 

RIN11930 

RIN11940 

RIN11950 

RIN11960 

RIN11970 

RIN11980 

RIN11990 

RIN12000 

RIN12010 

RIN12020 

RIN12030 

RIN12040 

RIN12050 

RIN12060 

RIN12070 

RIN12080 

■RIN12090 

RIN12100 

RIN12110 

RIN12120 

RIN12130 

RIN12140 

RIN12150 

RIN12160 

RIN12170 

RIN12180 

RIN12I90 

RIN12200 
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SUBROUTINE ELAS RIN12210 

C> — RIN12220 

C RIN12230 

C RIN12240 

C FUNCTION - DETERMINE ELASTIC SOLUTION RIN12250 

C (NO HYDRODYNAMICS) RIN12260 

C HELAS , HI ELAS , H2ELAS , H3ELAS , FORCE RIN12270 

C RIN12280 

C RESTRICTIONS - RIN12290 

C RIN12300 

C REMARKS - NOTE VARIABLES PASSED IN COMMON RIN12310 

RIN12320 

EXTERNAL REFERENCES - NONE RIN12330 

RIN12340 
RIN12350 

C> ; — — RIN12360 

IMPLICIT REAL*8 (A-H.O-Z) RIN12370 

DIMENSION YO(4) RIN12380 

COMMON/BPAR /AL, BET, PO ,EPS , PF RIN12390 

COMMON/BCD /XLOC( 100) ,C( 100) ,D( 100) ,NDX(3 ) ,NDP(3 ) ,DX(3 ) ,DP(3) RIN12400 

COMMON/BELAS /FORCE, HELAS, H1ELAS ,H2ELAS ,H3ELAS ,W1(4) ,W2(4) RIN12410 

COMMON/BFLAG /IFLAG, IPLOAD RIN12420 

YO(1)=0. RIN12430 

YO(2)=0. RIN12440 

YO(3)=0. RIN12450 

YO(4)=0. RIN12460 

WRITE(6 , 1009 ) RIN12470 

1009 F0RMAT(5X,1HX,8X,1HH,9X,2HH' ,7X,3HH' ' ,7X,4HH' ' ’ ,6X,4HPRES RIN12480 

$,8X,5HFORCE) RIN12490 

IF(IFLAG.EQ. 1 )THEN RIN12500 

IMIN=NDP(2)+1 RIN125 10 

IMAX=NDP(2)+NDP(3) RIN12520 . 

ELSE I F ( I FLAG . EQ . 2 ) THEN RIN12530 

IMIN=1 RIN12540 

IMAX=NDP(2) RIN12550 

ELSE IF(IFLAG.EQ.3)THEN RIN12560 

IMIN=1 RIN12570 

IMAX=NDP( 2 )+NDP( 3 ) RIN12580 

ELSE I F ( I FLAG . EQ . 4 ) THEN RIN12590 

IMIN=1 RIN12600 

IMAX=NDP ( 1 ) +NDP ( 2 ) RIN12610 

END IF RIN12620 

DO 600 I=IMIN, IMAX RIN12630 

YO( 1 )=YO( 1 )+C(I ) RIN12640 

YO(2)=YO(2)+D(I) RIN12650 

600 CONTINUE RIN12660 

YO(1)=YO(1)*PO*BET RIN12670 

YO( 2 )=YO( 2 )*P0*BET RIN12680 

HELAS=l.-YO(l) RIN12690 

H1ELAS=-Y0(2) RIN12700 

H2ELAS=-YO(3) RIN12710 

H3ELAS=-YO(4) RIN12720 

WRITE(6,'(F8.4,5F10.5,F13.5)')1., HELAS, H1ELAS,H2ELAS,H3ELAS,0. ,0. RIN12730 
FORCE=O.DO RIN12740 

IF(YO(l).GT.l.)THEN RIN12750 
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CHI=YO(1)-1.0 

Y0(l)=1.0 

B=CHI/ (W2(3)*W1( 1 )/Wl(3)-W2( 1 ) ) 

A=-B*W2(3)/W1(3) 

YO(2)=A*W1(2)+B*W2(2)+YO(2) 

Y0(3)=0. 

Y0(4)=A*W1(4)+B*W2(4) 

FORCE=YO ( 4 ) *AL/ BET 
HELAS=1.0-Y0(l) 

H1ELAS=-Y0(2) 

H2ELAS=-YO(3) 

H3ELAS--YO( 4 ) 

WRITE(6, '(F8.4,5F10.5,F13.5)')1.0,HELAS,H1ELAS,H2ELAS 
$ ,H3ELAS,0., FORCE 

END IF 
RETURN 
END 


RIN12760 

RIN12770 

RIN12780 

RIN12790 

RIN12800 

RIN12810 

RIN12820 

RIN12830 

RIN12840 

RIN12850 

RIN12860 

RIN12870 

RIN12880 

RIN12890 

RIN12900 

RIN12910 

RIN12920 


C>- 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c> 


SUBROUTINE ERRMSG 


FUNCTION - PRINT ERROR TERMS UPON EXIT FROM *ZSCNT* 

RESTRICTIONS - CALLED ONLY IF IER.NE.0 FROM *ZSCNT* 

REMARKS - FINAL VALUES OF X AND F FROM EVAL 

ARE PRINTED 

INPUT/OUTPUT: 

UNIT DESCRIPTION 

4 TERMINAL OUTPUT 

6 OUTPUT FILE 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/BLAST /X1,X2,F1,F2 

WRITE(4,'(8X,"W(1)" ,10X,"W(2)" ,10X,"F(D" ,10X,"F(2)' ' 

+ ,/)’) 

WRITE(4, ' (1X,4(2X,E12.5))' )X1,X2,F1,F2 

WRITE(6,'(8X,"W(D" ,10X,"W(2)' ',10X, ' ’ F( 1 ) ' ' , 10X, ' ’F(2)" 
+ )’) 

WRITE(6 , '(1X,4(2X,E12.5))')X1,X2,F1,F2 

RETURN 

END 


RIN12930 
RIN12940 
RIN12950 
RIN12960 
RIN12970 
RIN12980 
RIN12990 
RIN13000 
RIN13010 
RIN13020 
RIN13030 
RIN13040 
RIN13050 
RIN13060 
RIN13070 
RIN13080 
RIN13090 
RIN13100 
RIN131 10 
RIN13120 
RIN13130 
RIN13140 
RIN13150 
RIN13160 
RIN13170 
RIN13180 
RIN13190 
RIN13200 
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C> 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


c> 


5 


SUBROUTINE EVAL(X,F,N,PAR) 


FUNCTION - DEFINES EQUATIONS FOR H AND H’ WHICH ARE 

SOLVED BY SECANT METHOD *ZSCNT* 

RESTRICTIONS 

REMARKS - NOTE VARIABLE PASSED IN COMMON 

EXTERNAL REFERENCES - CONST 

INPUT/OUTPUT: 

UNIT DESCRIPTION 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

X X(l) IS DISPLACEMENT AT 1 

X(2 ) IS X' AT 1 

F EQUATIONS (2) FOR H AND H' 

N NO. OF EQUATIONS (N=2) 

PAR NOT USED 


IMPLICIT REAL*8 <A-H,0-Z) 

DIMENSION X(N),F(N),PAR(1) 

COMMON/BPAR /AL,BET,PO,EPS,PF 

COMMON/BCD /XLOC( 100) ,C( 100) ,D( 100) ,NDX(3 ) ,NDP(3 ) ,DX(3) , DP(3) 

COMMON/ BCOEFF /RK , RC , R1 , R2 , XCAV , ICAV 

COMMON/BINT /IELAS, IBACK, IPR 

COMMON/BLAST /X1LAST,X2LAST,F1LAST,F2LAST 

COMMON/BPROF /Cl 

COMMON/BFLAG /IFLAG, IPLOAD 

Hl=1.0-X(l) 

H11=C1-X(2) 

H0=H1-H11 

CALL CONST (HO, HI, PF) 

SUM1=X( 1 ) 

SUM2=X(2) 

INDEX=0 

CONTINUE 

INDEX=INDEX+1 

IF (I FLAG. EQ. 1 . AND. INDEX.GT.NDP(2 )+NDP(3 ) )GOTO 6 
IF ( I FLAG . EQ . 2 . AND . INDEX . GT . NDP ( 2 ) ) GOTO 6 
IF(XLOC( INDEX) . LT . DMIN1 ( 1 .-EPS , XCAV) )GOTO 5 
I F ( XLOC ( I NDEX ) . LT . XC AV ) TH EN 
PRES=0 
ELSE 

PRES=P(XLOC( INDEX)) 

END IF 

IF(XL0C(INDEX).GT.1.-EPS) PRES=PRES-PO 
CONTINUE 

SUM1=SUM1+C( INDEX )*PRES*BET 


RIN13210 

•RIN13220 

RIN13230 

RIN13240 

RIN13250 

RIN13260 

RIN13270 

RIN13280 

RIN13290 

RIN13300 

RIN13310 

RIN13320 

RIN13330 

RIN13340 

RIN13350 

RIN13360 

RIN13370 

RIN13380 

RIN13390 

RIN13400 

RIN13410 

RIN13420 

RIN13430 

RIN13440 

RIN13450 

RIN13460 

RIN13470 

RIN13480 

RIN13490 

RIN13500 

RIN13510 

RIN13520 

RIN13530 

RIN13540 

RIN13550 

RIN13560 

RIN13570 

RIN13580 

RIN13590 

RIN13600 

RIN13610 

RIN13620 

RIN13630 

RIN13640 

RIN13630 

RIN13660 

RIN13670 

RIN13680 

RIN13690 

RIN13700 

RIN13710 

RIN13720 

RIN13730 

RIN13740 

RIN13750 
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SUM2=SUM2+D( INDEX )*PRES*BET RIN1.3760 

GOTO 5 — RIN13770 

6 CONTINUE RIN13780 

F(1)=SUM1 RIN13790 

F(2)=SUM2 RIN13800 

X1LAST=X( 1 ) • RIN13810 

X2LAST=X(2 ) RIN13820 

F1LAST=F( 1 ) RIN13830 

F2LAST=F(2) RIN13840 

RETURN RIN13850 

END RIN13860 
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C> 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c> 


5 


SUBROUTINE EVAL3(X,F,N,PAR) 


FUNCTION - DEFINES EQUATIONS FOR H AND H' WHICH ARE 

SOLVED BY SECANT METHOD *ZSCNT* 

RESTRICTIONS 

REMARKS - NOTE VARIABLE PASSED IN COMMON 

EXTERNAL REFERENCES - CONST 

INPUT/OUTPUT: 

UNIT DESCRIPTION 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

X X(l) IS DISPLACEMENT AT 1 

X(2) IS X' AT 1 

F EQUATIONS (2) FOR H AND H' 

N NO. OF EQUATIONS (N=2) 

PAR NOT USED 


IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION X(N),F(N),PAR(1) 

COMMON/ BPAR /AL,BET,P0,EPS,PF 

COMMON/BCD /XLOC( 100) ,C( 100) , D( 100) ,NDX(3) ,NDP(3) ,DX(3) ,DP(3) 

COMMON/ BCOEFF / RK , RC , R 1 , R2 , XCAV , I CAV 

COMMON/BINT /IELAS , IBACK, IPR 

COMMON/BLAST /X1LAST,X2LAST,F1LAST,F2LAST 

COMMON/ BPROF / Cl 

COMMON/ BFLAG /IFLAG,IPL0AD 

Hl=1.0-X(l) 

H11=C1-X(2) 

H0=H1-H11 

CALL CONST ( HO, H1,PF) 

SUM1=X( 1 ) 

SUM2=X(2) 

INDEX=0 

CONTINUE 

INDEX=INDEX+1 

IF(INDEX.GT.NDP(2)+NDP(3))G0T0 6 
IF(XLOC(INDEX) .LT.DMIN1( 1 .-EPS ,XCAV ) )GOTO 5 
IF(XLOC(INDEX) .LT.XCAV)THEN 
PRES=0 
ELSE 

PRES=P(XLOC(INDEX) ) 

END IF 

PRES=PRES-PO 

CONTINUE 

SUM1=SUM1+C(INDEX)*PRES*BET 
SUM2=SUM2*D( INDEX )*PRES*BET 


RIN13870 
RIN13880 
RIN13890 
RIN13900 
RIN13910 
RIN13920 
RIN13930 
RIN13940 
RIN13950 
RIN13960 
RIN13970 
RIN13980 
RIN13990 
RIN14000 
RIN14010 
RIN14020 
RIN14030 
RIN14040 
RIN14050 
RIN14060 
RIN14070 
RIN14080 
RIN14090 
RIN14100 
RIN141 10 
RIN14120 
RIN14130 
RIN14140 
RIN14150 
RIN14160 
RIN14170 
RIN14180 
RIN14190 
RIN14200 
RIN14210 
RIN14220 
RIN14230 
RIN14240 
RIN14250 
RIN14260 
RIN14270 
RIN14280 
RIN14290 
RIN14300 
RIN14310 
RIN14320 
RIN14330 
RIN14340 
RIN14350 
RIN14360 
RIN14370 
RIN14380 
RIN14390 
RIN14400 
RIN14410 
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GOTO 5 

RIN14420 

CONTINUE 

RIN14430 

F( 1 )=SUM1 

RIN14440 

F(2)=SUM2 

RIN14450 

X1LAST=X( 1) 

RIN14460 

X2LAST=X(2) 

RIN14470 

F1LAST=F( 1 ) 

RIN14480 

F2LAST=F(2) 

RIN14490 

RETURN 

RIN14500 

END 

RIN14510 
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SUBROUTINE EVAL4(X,F,N,PAR) 

C> 

C 

c 

C FUNCTION - DEFINES EQUATIONS FOR H AND H' WHICH ARE 

C SOLVED BY SECANT METHOD *ZSCNT* 

C 

C RESTRICTIONS 

C 

C REMARKS - NOTE VARIABLE PASSED IN COMMON 

C 

C EXTERNAL REFERENCES - CONST 

C 

C INPUT/OUTPUT: 

C UNIT DESCRIPTION 

C 

C ARGUMENT DEFINITION: 

C NAME DESCRIPTION 

C X X(l) IS DISPLACEMENT AT 1 

C X(2) IS X' AT 1 

C F EQUATIONS (2) FOR H AND H* 

C N NO. OF EQUATIONS (N=2) 

C PAR NOT USED 

C 

C 

c > 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION X(N),F(N),PAR(1) 

COMMON/ BPAR / AL , BET , PO , EPS , PF 

COMMON/BCD /XLOC( 100) ,C( 100) ,D(100) ,NDX(3) ,NDP(3) ,DX(3) ,DP(3) 

COMMON/ BCOEFF /RK, RC, R1 ,R2 ,XCAV, I CAV 

COMMON/ BINT /IELAS,IBACK,IPR 

COMMON/ BLAST /X1LAST,X2LAST,F1LAST,F2LAST 

COMMON/ BPROF /Cl 

COMMON/ BFLAG /I FLAG, IPLOAD 

Hl=1.0-X(l) 

H11=C1-X(2) 

H0=H1-H11 

CALL CONST(HO ,H1 ,PF ) 

SUM1=X( 1 ) 

SUM2=X(2) 

INDEX=0 

5 CONTINUE 

INDEX=INDEX+1 

IF( INDEX. GT.NDP( 1 )+NDP(2) )GOTO 6 

IF(XLOC( INDEX ) ,LT.DMIN1( 1 . -EPS, XCAV))GOTO 5 ' 

I F ( XLOC ( INDEX ) . LT . XCAV )THEN 
PRES=0 
ELSE 

PRES=P(XLOC( INDEX)) 

END IF 

PRES=PRES-PO 
7 CONTINUE 

SUM1=SUM1+C(INDEX)*PRES*BET 
SUM2=SUM2+D( INDEX)*PRES*BET 


RIN14520 

RIN14530 

RIN14540 

RIN14550 

RIN14560 

RIN14570 

RIN14580 

RIN14590 

RIN14600 

RIN14610 

RIN14620 

RIN14630 

RIN14640 

RIN14650 

RIN14660 

RIN14670 

RIN14680 

RIN14690. 

RIN14700 

RIN14710 

RIN14720 

RIN14730 

RIN14740 

RIN14750 

RIN14760 

RIN14770 

RIN14780 

RIN14790 

RIN14800 

RIN14810 

RIN14820 

RIN14830 

RIN14840 

RIN14850 

RIN14860 

RIN14870 

RIN14880 

RIN14890 

RIN14900 

RIN14910 

RIN14920 

RIN14930 

RIN14940 

RIN14950 

RIN14960 

RIN14970 

RIN14980 

RIN14990 

RIN15000 

RIN15010 

RIN15020 

RIN15030 

RIN15040 

RIN15050 

RIN15060 
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GOTO 5 

6 CONTINUE 
F ( 1 )=SUM1 
F(2)=SUM2 
XlLAST=X( 1 ) 
X2LAST=X(2) 
F1LAST=F( 1 ) 
F2LAST=F(2) 
RETURN 
END 


RIN15070 

RIN15080 

RIN15090 

RIN15100 

RIN15110 

RIN15120 

RIN15130 

RIN15140 

RIN15150 

RIN15160 


C> 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C> 


FUNCTION P(PSI) 


FUNCTION - RETURN HYDRODYNAMIC PRESSURE AT PSI 

RESTRICTIONS 

REMARKS - NOTE VARIABLES PASSED IN COMMON 

EXTERNAL REFERENCES - DABS 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

PSI VALUE OF X WHERE PRESSURE IS DESIRED 


IMPLICIT REAL*8 (A-H.O-Z) 

COMMON / BCOEFF / RK , RC , R 1 , R2 , XCAV , I C AV 

COMMON/BINT /IELAS,IBACK,IPR 

H=R1+R2*PSI 

HINV=1./H 

P=0. 

IF(PSI. GE. XCAV. AND. DABS(R2).GT.1.D-10)THEN 
P=1 . /R2*(HINV*(-1 .+RK/2 .*HINV)+RC ) 

END IF 

I F ( I BACK . EQ . 1 ) THEN 
P=-P 
END IF 
RETURN 
END 


RIN15170 

-RIN15180 

RIN15190 

RIN15200 

RIN15210 

RIN15220 

RINI5230 

RIN15240 

RIN15250 

RIN15260 

RIN15270 

RIN15280 

RIN15290 

RIN15300 

RIN153IO 

RIN15320 

RIN15330 

-RIN15340 

RIN15350 

RIN15360 

RIN15370 

RIN15380 

RIN15390 

RIN15400 

RIN15410 

RIN15420 

RIN15430 

RIN15440 

RINI5450 

RIN15460 

RINI5470 

RINI5480 
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SUBROUTINE PRT(X,Y, ITERM) 

C> 

c 

c 

C FUNCTION - USING DISPLACEMENTS, CONVERT TO FILM 

C THICKNESSES AND PRINT OUT 

C 

C RESTRICTIONS 

C 

C REMARKS 

C 

C EXTERNAL REFERENCES - P 

C 

C INPUT/OUTPUT: 

C UNIT DESCRIPTION 

C ITERM OUTPUT UNIT 

C 

C ARGUMENT DEFINITION: 

C NAME DESCRIPTION 

C X DIMENSIONLESS LENGTH VARIABLE . 

C Y DISPLACEMENT AND ITS DERIVATIVES 

C ITERM OUTPUT UNIT 

C 
C 

C> -- — 

IMPLICIT REAL*8 (A-H.O-Z) , ■ V 

COMMON/BINT /IELAS,IBACK,IPR 
COMMON/BPROF /Cl 
DIMENSION Y(4) 

H=1.-Y(1)+C1*(X-1.) 

H1=C1-Y(2) ' , . ' ... 

H2=-Y(3) . : . ‘ 

H3=-Y(4) ! 

PRES=0 . 

IF(IELAS.EQ.O.AND.X.GE.O. )PRES=P(X) 
WRITE(ITERM,5)X,H,H1,H2,H3, PRES 
5 F0RMAT(F10.4,5F12.5) 

RETURN 

END 


RIN15490 

-RIN15500 

RIN15510 

RIN15520 

RIN15530 

RIN15540 

RIN15550 

RIN15560 

RIN15570 

RIN15580 

RIN15590 

RIN15600 

RIN15610 

RIN15620 

RIN15630 

RIN15640 

RIN15650 

RIN15660 

RIN15670 

RIN15680 

RIN15690 

RIN15700 

RIN15710 

RIN15720 

-RIN15730 

RIN15740 

RIN15750 

RIN15760 

RIN15770 

RIN15780 

RIN15790 

RIN15800 

RIN15810 

RIN15820 

RIN15830 

RIN15840 

RIN15850 

RIN15860 

RIN15870 
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C> 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c> 


c 


SUBROUTINE PRTOUT(XX) 


FUNCTION - PRINT FILM THICK., PRES., ETC. 

RESTRICTIONS 

REMARKS - NOTE VARIABLES PASSED IN COMMON 

EXTERNAL REFERENCES - DFN I , DFN2 , DFN3 , PRT , RUK 

INPUT/OUTPUT: 

UNIT DESCRIPTION 

6 OUTPUT FILE 

ARGUMENT DEFINITION: 

NAME DESCRIPTION 

XX(2) DISPLACEMENT, SLOPE OF DISPLACEMENT 


IMPLICIT REAL*8 (A-H,0-Z) 

COMMON/BINT /IELAS , IBACK, IPR 

COMMON/BELAS /FORCE , HELAS , HlELAS , H2ELAS , H3ELAS , W1 ( 4 ) , W2 ( 4 ) 
COMMON/BCD /XLOC(100),C(100),D(100),NDX(3),NDP(3),DX(3),DP(3) 
COMMON / B FLAG /IFLAG, IPLOAD 

DIMENSION XX(2) , YO(4), YT(4) ,DJ(4),CKJ(4,4) 

EXTERNAL DFN1 ,DFN2 ,DFN3 


XNN=1 .0 
YO( 1 )=XX( 1 ) 

YO(2)=XX(2) 

Y0(3)=0. 

YO(4)=0. 

I F ( IBACK . EQ . 1 . AND . FORCE . GT . I . D-8 )THEN 
WRITE(6 , 1009 ) 

WRITE(6, '(F8.4,5F10.5,F13.5)' )1.0, HELAS, HlELAS, H2ELAS 
$ ,H3ELAS,0. , FORCE 

RETURN 
ELSE 

WRITE(6 , 1007 ) 

IF(IPR.EQ. 1 )THEN 

WRITE(6 , ' ( ) 

ELSE 

MRITE(6,'(” ”)’) 

END IF 

CALL PRT(XNN,YO,6) 

END IF 

IF(IPR.EQ. 1 )THEN 

IF( IFLAG. EQ. 2. OR. IFLAG. EQ.4)G0T0 111 
IPLOAD=l 

CALL RUK ( -DX ( 3 ) / 2 . , -DP( 3 ) /2 . , XNN , XNN , YO , YO , 4 , DFN2 , YT , DJ , CK J , 
CALL PRT(XNN,YO,6) 

DO 9 JJ=2 ,NDP(3 ) 


RIN15880 

■ RIN15890 

RIN15900 

RIN15910 

RIN15920 

RIN15930 

RIN15940 

RIN15950 

RIN15960 

RIN15970 

RIN15980 

RIN15990 

RIN16000 

RIN16010 

RIN16020 

RIN16030 

RIN16040 

RIN16050 

RIN16060 

RIN16070 

RIN16080 

RIN16090 

RIN16100 
RIN161 10 
RIN16120 
RIN16130 
RIN16140 
RIN16150 
RIN16160 
RIN16170 
RIN16180 
RIN16190 
RIN16200 
RIN16210 
RIN16220 
RIN16230 
RIN16240 
RIN16250 
RIN16260 
RIN16270 
RIN16280 
RIN16290 
RIN16300 
RIN16310 
RIN16320 
RIN16330 
RIN16340 
RIN16350 
RIN16360 
RIN16370 
RIN16380 
RIN16390 
4(), RIN16400 
RIN16410 
RIN16420 
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J=NDP(3)+1-JJ 

CALL RUK(-DX(3),-DP(3),XNN,XNN,YO,YO,4,DFN2,YT,DJ,CKJ,4) 

CALL PRT(XNN,YO,6) 

CONTINUE 

CALL RUK(-DX(3)/2.,-DP(3)/2.,XNN,XNN,YO,YO,4,DFN2,YT,DJ,CKJ,4) 
CALL PRT(XNN,YO,6) 

CONTINUE 

IF(IFLAG.EQ.1)THEN 

IPLOAD=0 

ELSE IF(IFLAG.EQ.2)THEN 
IPLOAD=l 

ELSE IF(IFLAG.EQ.3)THEN 
IPLOAD=l 

ELSE IF(IFLAG.EQ.4)THEN 
IPLOAD=l 
END IF 

CALL RUK(-DX(2)/2.,-DP(2)/2.,XNN,XNN,YO,YO,4,DFN2,YT,DJ,CKJ,4) 
CALL PRT(XNN,YO,6) 

DO 8 JJ=2 ,NDP(2 ) 

J=NDP(2 )+l-JJ 

CALL RUK( -DX ( 2 ) , -DP { 2 ) , XNN , XNN , YO , YO , 4 , DFN2 , YT , D J , CK J , 4 ) 

CALL PRT(XNN,YO,6) 

CONTINUE 

CALL RUK(-DX(2)/2.,-DP(2)/2.,XNN,XNN,YO,YO,4,DFN2,YT,DJ,CKJ > 4) 
CALL PRT ( XNN , YO , 6 ) 

IF(IFLAG.EQ. 1.0R.IFLAG.EQ.2)THEN 
WRITE(6, 1 (' '(§")' ) 

RETURN 
END IF 

DO 7 J=1,NDP(1) 

IF(IFLAG.EQ.4)THEN 

CALL RUK (-DX(1),-DP(1), XNN , XNN , YO , YO , 4 , DFN3 , YT , D J , CK J , 4 ) 
ELSE 

CALL RUK(-DXd),-DPd),XNN,XNN,YO,YO,4,DFNl,YT,DJ,CKJ,4) 

END IF 

CALL PRT(XNN,YO,6) 

CONTINUE 
END IF 

WRITE(6, ' ( "Q 1 1 ) ' ) 

RETURN 


1007 FORMAT(7X,1HX,10X,1HH,11X,2HH' ,9X,3HH' ' ,9X,4HH' ' ' ,8X,4HPRES) 
1009 F0RMAT(5X,1HX,8X,1HH,9X,2HH' ,7X,3HH' ' ,7X,4HH’ ' ' ,6X,4HPRES 
$,8X,5HFORCE) 

C 

END 
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RIN16890 




c 


RIN16910 

c 

• 

RIN16920 

c 

FUNCTION - RUNGE-KUTTA INTEGRATION 

RIN16930 

c 


- RIN16940 

c 

RESTRICTIONS 

RIN16950 

c 


RIN16960 

c 

REMARKS 

RIN16970 

c 


RIN16980 

c 

EXTERNAL REFERENCES - DFN 

RIN16990 

c 


RIN17000 

c 

ARGUMENT DEFINITION: 

RIN17010 

c 

NAME DESCRIPTION 

RIN17020 

c 

1 ' • • - * . . 

RIN17030 

c 


RIN17040 

DTX11 7H Q n 





IMPLICIT REAL* 8 (A-H,0-Z) 

RIN17060 


DIMENSION YO(ID) , YN(ID) ,D(ID) ,YT(ID) ,CK(4,ID) ' 

RIN17070 


N=DP/DX+.001 

RIN17080 


XN=XO 

RIN17090 


DO 2 J=l,NO 

RIN17100 


2 YN( J)=YO( J) 

RIN17110 


DO 5 1=1, N 

RIN17120 


CALL DFN(XN,YN,D) 

RIN17130 


DO 6 J=1 ,NO 

RIN17140 


CK(1,J)=D(J)*DX 

RIN17150 


6 YT(J)=YN( J)+CK( 1 , J)/2. 

RIN17160 


CALL DFN(XN+DX/2. ,YT,D) 

RIN17170 


DO 7 J=l,NO 

RIN17180 


CK(2,J)=D(J)*DX 

RIN17190 


7 YT(J)=YN(J)+CK(2,J)/2. ‘ 

RIN17200 


CALL DFN ( XN+DX / 2 . , YT , D ) ' 

RIN17210 


DO 8 J=l,NO 

RIN17220 


CK(3,J)=D(J)*DX 

RIN17230 


8 YT( J)=YN( J)+CK(3, J) 

RIN17240 


CALL DFN ( XN+DX , YT , D ) 

RIN17250 


DO 9 J=l,NO 

RIN17260 


CK(4, J)=D(J)*DX . 

RIN17270 


YN(J)=YN(J)+(CK(l,J)+2.*(CK(2,J)+CK(3,J))+CK(4,J))/6. 

RIN17280 

9 

CONTINUE 

RTN17290 


5 XN=XN+DX , 

RIN17300 


RETURN 

RIN17310 


END. 

RIN17320 
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C 

C 

C- 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE ZSCNT (FCN, NSIG.N, ITMAX, PAR, X, FNORM, WK, IER) 
IMSL ROUTINE NAME - ZSCNT 

THIS PROGRAM SHOULD BE COMPILED IN FORTRAN 4 (FORTHX) 


COMPUTER - IBM/DOUBLE 

LATEST REVISION - JUNE 1, 1980 

PURPOSE - SOLVE A SYSTEM OF NONLINEAR EQUATIONS 


USAGE 


- CALL ZSCNT (FCN, NSIG.N, ITMAX,PAR,X,FNORM, 
WK.IER) 


ARGUMENTS FCN - THE NAME OF A USER-SUPPLIED SUBROUTINE WHICH 

EVALUATES THE SYSTEM OF EQUATIONS TO BE 
SOLVED. FCN MUST BE DECLARED EXTERNAL IN 
THE CALLING PROGRAM AND MUST HAVE THE 
FOLLOWING FORM, 

SUBROUTINE FCN(X,F,N,PAR) 

DIMENSION X(N) ,F(N) ,PAR(1) 

F(l)= 


F(N)= 

RETURN 

END 

GIVEN X(1)...X(N), FCN MUST EVALUATE THE 
FUNCTIONS F(1)...F(N) WHICH ARE TO BE MADE 
ZERO. X SHOULD NOT BE ALTERED BY FCN. THE 
PARAMETERS IN VECTOR PAR (SEE ARGUMENT 
PAR BELOW) MAY ALSO BE USED IN THE 
CALCULATION OF F(1)...F(N). 

NSIG - THE NUMBER OF DIGITS OF ACCURACY DESIRED 
IN THE COMPUTED ROOT (INPUT). 

N - THE NUMBER OF EQUATIONS TO BE SOLVED AND 

THE NUMBER OF UNKNOWNS (INPUT). 

ITMAX - THE MAXIMUM ALLOWABLE NUMBER OF ITERATIONS 
(INPUT). 

PAR - PAR CONTAINS A PARAMETER SET WHICH IS 

PASSED TO THE USER SUPPLIED FUNCTION FCN. 
PAR MAY BE USED TO PASS ANY AUXILIARY 
PARAMETERS NECESSARY FOR COMPUTATION OF 
THE FUNCTION FCN. (INPUT) 

X - A VECTOR OF LENGTH N. (INPUT/ OUTPUT) ON INPUT 

X IS THE INITIAL APPROXIMATION TO THE ROOT. 
ON OUTPUT, X IS THE BEST APPROXIMATION TO 
THE ROOT FOUND BY ZSCNT. 

FNORM - ON OUTPUT, FNORM IS EQUAL TO 

F(l)**2+. . .F(N)**2 AT THE POINT X. 

WK - WORK VECTOR OF LENGTH (N+l)*(3*N+8) 

IER - ERROR PARAMETER. (OUTPUT) 

TERMINAL ERROR 

IER = 129 INDICATES THAT ZSCNT FAILED TO 
CONVERGE WITHIN ITMAX ITERATIONS. THE 
USER MAY INCREASE ITMAX OR TRY A NEW 
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C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

G 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 

c 

c 


PRECISION/HARDWARE 
REQD. IMSL ROUTINES 


NOTATION 

COPYRIGHT 

WARRANTY 


INITIAL GUESS. 

IER = 130 INDICATES THE ALGORITHM WAS 
UNABLE TO IMPROVE ON THE RETURNED VALUE 
OF X. THIS SITUATION ARISES WHEN THE 
SOLUTION CANNOT BE DETERMINED TO NSIG 
DIGITS DUE. TO ERRORS IN THE FUNCTION 
VALUES. IT MAY ALSO INDICATE THAT THE 
ROUTINE IS TRAPPED IN THE AREA OF A 
LOCAL MINIMUM. THE USER MAY TRY A NEW 
INITIAL GUESS. 

- SINGLE AND DOUBLE/H32 

- SINGLE/H36 ,H48 ,H60 

- S INGLE/ GGUBFS , LEQT2F , LUDATF , LUELMF , LUREFF , 

UERSET , UERTST , UGETIO , ZSCNU 

- DOUBLE /GGUBFS , LEQT2F , LUDATF , LUELMF , LUREFF , 

UERSET , UERTST , UGETIO , VXADD , VXMUL , VXSTO , 

ZSCNU 

- INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 

INTRODUCTION OR THROUGH IMSL ROUTINE UHELP 

- 1980 BY IMSL, INC. ALL RIGHTS RESERVED. 

- IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 

EXPRESSED OR IMPLIED, IS APPLICABLE. 


INTEGER 

DOUBLE PRECISION 

INTEGER 

EXTERNAL 


SUBROUTINE ZSCNT (FCN, NSIG, N,ITMAX, PAR, X,FNORM,WK, IER) ; 

SPECIFICATIONS FOR ARGUMENTS 
IER, ITMAX,N, NSIG 
FNORM, PARC 1) ,WK(1) ,X(N) 

SPECIFICATIONS FOR LOCAL VARIABLES 
11,12,13,14, 15 ,LNEW,L0LD,N1 
FCN 

FIRST EXECUTABLE STATEMENT 

N1=N+1 

11 = N1*N1+1 

12 = I1+N*N1 

13 = I2+N1 

14 = I3+N1 

15 = I4+-N1 

CALL UERSET (O.LOLD) 

CALL ZSCNU(X,N, FCN, NSIG, N1,WK(1) ,WK(I1) ,WK(I2) ,WK(I3) ,WK(I4) 

* ,WK(I5) , ITMAX,PAR, IER) 

CALL FCN(X,WK,N,PAR) 

FNORM = 0.0D0 
DO 5 I = 1,N 

FNORM = FNORM + WK(I)*WK(I) 

5 CONTINUE 

CALL UERSET (LOLD,LNEW) 
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IF (IER.EQ.O) GO TO 9005 

ZSC01110 

9000 

CONTINUE 

ZSC01120 


CALL UERTST ( IER , 6HZSCNT ) 

ZSC01130 

9005 

RETURN 

ZSC01140 


END 

ZSC01150 


FUNCTION GGUBFS (DSEED) 

ZSC01160 


DOUBLE PRECISION DSEED 

ZSC01170 


DOUBLE PRECISION D2P31M.D2P31 

ZSC01180 


DATA D2P31M/2147483647 .DO/ 

ZSC01190 


DATA D2P31 /214 7483648 . DO/ 

ZSC01200 


DSEED = DMOD ( 16807 . DO*DSEED ,D2P3 1M) 

ZSC01210 


GGUBFS = DSEED / D2P31 

ZSC01220 


RETURN 

ZSC01230 


END 

ZSC01240 
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SUBROUTINE LEQT2F (A, M, N , IA,B , IDGT .WKAREA, IER) 

DIMENSION A(IA,1),B(IA,1),WKAREA(1) 

DOUBLE PRECISION A, B, WKAREA, D1.D2.WA 

IER=0 

JER=0 

J = N*N+1 

K = J+N 

MM = K+N 

KK = 0 

MM1 = MM-1 

JJ=1 

DO 5 L=1,N 
DO 5 1=1, N 

WKAREA (JJ)=A( I, L) 

JJ=JJ+1 
5 CONTINUE 

CALL LUDATF (WKAREA(l) , A,N ,N , IDGT ,D1 ,D2 ,WKAREA( J) .WKAREA(K) , 

1 WA.IER) 

IF (IER.GT. 128) GO TO 25 

IF (IDGT .EQ. 0 .OR. IER .NE. 0) KK = 1 

DO 15 I = 1,M 

CAllL LUELMF (A,B(1,I) .WKAREA (J) ,N ,N ,WKAREA(MM) ) 

IF (KK .NE. 0) 

* CALL LUREFF (WKAREA(l) ,B(1, I) ,A,WKAREA(J) ,N,N, WKAREA (MM) ,IDGT, 

* WKAREA(K) ,WKAREA(K) , JER) 

DO 10 11=1, N 

B ( I I , I) = WKAREA (MM1+I I) 

10 CONTINUE 

IF (JER.NE.O) GO TO 20 
15 CONTINUE 
GO TO 25 
20 IER = 131 
25 JJ=1 

DO 30 J = 1 ,N 
DO 30 I = 1,N 

A(I , J)=WKAREA( JJ) 

JJ=JJ+1 
30 CONTINUE - 

IF (IER .EQ. 0) GO TO 9005 
9000 CONTINUE 

CALL UERTST (IER.6HLEQT2F) 

9005 RETURN 
' END 
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ZSC01340 
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SUBROUTINE LUDATF (A,LU,N, IA, IDGT.Dl ,D2, IPVT,EQUIL,WA, 

IER) 

ZSC01680 

DIMENSION A(IA, 1) ,LU( IA, 1) , IPVT(l) ,EQUIL(1) 


ZSC01690 

DOUBLE PRECISION A,LU,D1 ,D2 .EQUIL.WA, ZERO, ONE, FOUR, SIXTN.SIXTH, 

ZSC01700 

* 

RN, WREL, BIGA, BIG, P, SUM, AI.WI.T, TEST, Q 

ZSC01710 

DATA ZERO, ONE, FOUR, SIXTN, SIXTH/O. DO, 1. DO 

,4. DO, 

ZSC01720 

* 

16. DO, . 0625D0/ 


ZSC01730 

IER = 0 


ZSC01740 

RN 

= N 


ZSC01750 

WREL = ZERO 


ZSC01760 

D1 

= ONE 


ZSC01770 

D2 

= ZERO 


ZSC01780 

BIGA = ZERO 


ZSC01790 

DO 

10 1=1, N 


ZSC01800 


BIG = ZERO 


ZSC01810 


DO 5 J=1,N 


ZSC01820 


P = A(I,J) 


ZSC01830 


LU(I , J) = P 


ZSC01840 


P = DABS (P) 


ZSC01850 


IF (P .GT. BIG) BIG = P 


ZSC01860 

5 

CONTINUE 


ZSC01870 


IF (BIG .GT. BIGA) BIGA = BIG 


ZSC01880 


IF (BIG .EQ. ZERO) GO TO 110 


ZSC01890 


EQUIL(I) = ONE/BIG 


ZSC01900 

10 CONTINUE 


ZSC019 10 

DO 

105 J=1,N 


ZSC01920 


JM1 = J-l 


ZSC01930 


IF ( JM1 . LT. 1) GO TO 40 


ZSC01940 


DO 35 1=1, JM1 


ZSC01950 


SUM = LU(I , J) 


ZSC01960 


IM1 = 1-1 


ZSC01970 


IF (IDGT .EQ. 0) GO TO 25 


ZSC01980 


AI = DABS (SUM) 


ZSC01990 


WI = ZERO 


ZSC02000 


IF (IM1 .LT. 1) GO TO 20 


ZSC02010 


DO 15 K=1 , IM1 


ZSC02020 


T = • LU(I ,K)*LU(K, J) 


ZSC02030 


SUM = SUM-T 


ZSC02040 


WI = WI+DABS(T) 


ZSC02050 

15 

CONTINUE 


ZSC02060 


LU(I , J) = SUM 


ZSC02070 

20 

WI = WI+DABS(SUM) 


ZSC02080 


IF (AI .EQ. ZERO) AI = BIGA 


ZSC02090 


TEST = WI/AI 


ZSC02100 


IF (TEST .GT. WREL) WREL = TEST 


ZSC02110 


GO TO 35 


ZSC02120 

25 

IF ( IM1 .LT. 1) GO TO 35 


ZSC02130 


DO 30 K=1 , IM1 


ZSC02140 


SUM = SUM-LU(I ,K)*LU(K, J) 


ZSC02150 

30 

CONTINUE 


ZSC02160 


LU(I , J) = SUM 


ZSC02170 

35 

CONTINUE 


ZSC02180 

40 

P = ZERO 


ZSC02190 


DO 70 I=J,N 


ZSC02200 


SUM = LU(I , J) 


ZSC02210 


IF (IDGT .EQ. 0) GO TO 55 


ZSC02220 
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AI = DABS (SUM) 

ZSC02230 


WI = ZERO 

ZSC02240 


IF ( JM1 .LT. 1) GO TO 50 

ZSC02250 


DO 45 K= 1 , JM 1 

ZSC02260 


T = LU(I,K)*LU(K,J) 

ZSC02270 


SUM = SUM-T 

ZSC02280 


WI = WI+DABS (T) 

ZSC02290 

45 

CONTINUE 

ZSC02300 


LU(I , J) = SUM 

ZSC02310 

50 

WI = WI+DABS (SUM) 

ZSC02320 


IF (AI .EQ. ZERO) AI = BIGA 

ZSC02330 


TEST = WI/AI 

ZSC02340 


IF (TEST .GT. WREL) WREL = TEST 

. ZSC02350 


GO TO 65 

ZSC02360 

55 

IF (JM1 .LT. 1) GO TO 65 

ZSC02370 


DO 60 K= 1 , JM 1 

ZSC02380 


SUM = SUM-LU(I ,K)*LU(K, J) 

ZSC02390 

60 

CONTINUE 

ZSC02400 


LU(I , J) = SUM 

ZSC02410 

65 

Q = EQUIL(I)*DABS (SUM) 

ZSC02420 


IF (P .GE. Q) GO TO 70 

ZSC02430 


P = Q 

ZSC02440 


IMAX = I 

ZSC02450 

70 

CONTINUE 

ZSC02460 


IF (RN+P .EQ. RN) GO TO 110 

ZSC02470 


IF (J .EQ. IMAX) GO TO 80 

ZSC02480 


D1 = -D1 

ZSC02490 


DO 75 K=1,N 

ZSC02500 


P = LU(IMAX,K) 

ZSC02510 


LU(IMAX,K) = LU(J,K) 

ZSC02520 


LU( J,K) = P 

ZSC02530 

75 

CONTINUE 

ZSC02540 


EQUIL(IMAX) = EQUIL(J) 

ZSC02550 

80 

IPVT(J) = IMAX 

ZSC02560 


D1 = D1*LU(J, J) 

ZSC02570 

85 

IF (DABS (Dl) .LE. ONE) GO TO 90 

ZSC02580 


D1 = D1*SIXTH 

ZSC02590 


D2 = D2+F0UR 

ZSC02600 


GO TO 85 

ZSC02610 

90 

IF (DABS (Dl) .GE. SIXTH) GO TO 95 

ZSC02620 


Dl = D1*SIXTN 

ZSC02630 


D2 = D2-FQUR 

ZSC02640 


GO TO 90 

ZSC02650 

95 

CONTINUE 

ZSC02660 


JP1 = J+l 

ZSC02670 


IF (JP1 .GT. N) GO TO 105 

ZSC02680 


P = LU(J, J) 

ZSC02690 


DO 100 I=JP1 ,N 

ZSC02700 


LU(I , J) = LU(I,J)/P 

ZSC027 10 

100 

CONTINUE 

ZSC02720 

105 

CONTINUE 

ZSC02730 


IF (IDGT .EQ. 0) GO TO 9005 

ZSC02740 


P = 3*N+3 

ZSC02750 


WA = P*WREL 

ZSC02760 


IF (WA+10 .D0**(-IDGT) .NE. WA) GO TO 9005 

ZSC027 70 
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IER = 34 
GO TO 9000 
110 IER = 129 
D1 = ZERO 
D2 = ZERO 
9000 CONTINUE 

CALL UERTST ( IER , 6HLUDATF ) 
9005 RETURN 
END 


ZSC02780 

ZSC02790 
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ZSC02810 

ZSC02820 

ZSC02830 

ZSC02840 

ZSC02850 

ZSC02860 


SUBROUTINE LUELMF (A,B , IPVT.N, IA,X) 

DIMENSION A(IA,1),B(1),IPVT(1),X(1) 

DOUBLE PRECISION A,B,X,SUM 
DO 5 1=1, N 
5 X(I) = B(I) 

IW = 0 
DO 20 1=1, N 
IP = IPVT(I) 

SUM = X(IP) 

X(IP) = X(I) 

IF (IW .EQ. 0) GO TO 15 
IM1 = 1-1 
DO 10 J=IW, IM1 

SUM = SUM-A(I,J)*X(J) 

10 CONTINUE 

GO TO 20 

15 IF (SUM .NE. 0 .DO) IW = I 

20 X(I) = SUM 
DO 30 IB=1,N 
I = N+l-IB 
IP1 = 1+1 
SUM = X(I) 

IF (IP1 .GT. N) GO TO 30 
DO 25 J=IP1,N 

SUM = SUM-A(I , J)*X(J) 

25 CONTINUE 
30 X(I) = SUM/A(I , I) 

RETURN 

END 
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SUBROUTINE LUREFF (A,B ,UL, IPVT.N, IA,X, IDGT ,RES ,DX, IER) 


A(IA,1),UL(IA,1),B(1),X(1) ,RES(1),DX(1), 
ACCXTC2) 

A , ACCXT , B , UL , X , RES , DX , ZERO , XNQRM , DXNORM 
ITMAX/ 75/, ZERO/O. DO/ 


DIMENSION 
DIMENSION 
DOUBLE PRECISION 
DATA 
IER=0 

XNORM = ZERO 
DO 10 1=1, N. 

XNORM = DMAX1 (XNORM, DABS (X(I))) 

10 CONTINUE 

IF (XNORM .NE. ZERO) GO TO 20 
IDGT = 50 
GO TO 9005 

20 DO 45 ITER*1, ITMAX 
DO 30 1=1, N 
ACCXT ( 1 ) = 0.0D0 
ACCXT(2) = 0.0D0 

CALL VXADD(B(I) , ACCXT) 

DO 25 J=1,N 

CALL VXMUL(-A(I,J) ,X(J) .ACCXT) 

25 CONTINUE 

CALL VXSTO(ACCXT,RES(I) ) 

30 CONTINUE 

CALL LUELMF (UL, RES , IPVT.N, IA.DX) 

DXNORM = ZERO 
XNORM = ZERO 
DO 35 1=1, N 

X(I) = X(I) + DX(I) 

DXNORM = DMAX1 (DXNORM , DABS (DX ( I ) ) ) 

XNORM = DMAX1 (XNORM, DABS (X( I) ) ) 

35 CONTINUE 

IF (ITER .NE. 1) GO TO 40 
IDGT =50 

IF (DXNORM .NE. ZERO) IDGT = - DLOG 1 0 ( DXNORM / XNORM ) 
40 IF (XNORM+DXNORM .EQ. XNORM) GO TO 9005 

45 CONTINUE 
IER = 129 
9000 CONTINUE 

CALL UERTST ( IER , 6HLUREFF ) 

9005 RETURN 
END 
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1 SUBROUTINE UERSET (LEVEL, LEVOLD) 
INTEGER LEVEL, LEVOLD 

LEVOLD = LEVEL 

CALL UERTST (LEVOLD, 6HUERSET) 

RETURN 

END 


ZSC03570 

ZSC03580 

ZSC03590 

ZSC03600 

ZSC03610 
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SUBROUTINE UERTST (IER.NAME) 

INTEGER IEQ 

INTEGER* 2 - , NAME ( 3 ) , NAMSET ( 3 ) , NAMEQ ( 3 ) 

DATA NAMSET/ 2HUE , 2HRS , 2HET/ 

DATA NAMEQ/ 2H ,2H ,2H / 

DATA LEVEL/4/, IEQDF/O/, IEQ/ 1H=/ 

IF (IER.GT. 999) GO TO 25 

IF (IER.LT.-32) GO TO 55 

IF (IER.LE. 128) GO TO 5 

IF (LEVEL.LT. 1) GO TO 30 

CALL UGETI0(1 ,NIN, IOUNIT) 

IF (IEQDF.EQ. 1) WRITE (IOUNIT ,35) IER, NAMEQ, IEQ, NAME 
IF (IEQDF.EQ. 0) WRITE (IOUNIT, 35) IER.NAME 
GO TO 30 

5 IF (IER.LE;. 64) GO TO 10 
IF (LEVEL.LT. 2) GO TO 30 
CALL UGETI0(1,NIN, IOUNIT) 

IF (IEQDF.EQ. 1) WRITE (IOUNIT, 40) IER, NAMEQ, IEQ, NAME 
IF (IEQDF.EQ. 0) WRITE (IOUNIT, 40) IER.NAME 
GO TO 30 

10 IF (IER.LE. 32) GO TO 15 
IF (LEVEL.LT. 3) GO TO 30 
CALL UGETI0(1,NIN, IOUNIT) 

IF (IEQDF.EQ. 1) WRITE (IOUNIT, 45) IER, NAMEQ, IEQ, NAME 
IF (IEQDF.EQ. 0) WRITE (IOUNIT, 45) IER.NAME 
GO TO 30 
15 CONTINUE 
DO 20 1=1,3 

IF (NAME ( I ) . NE . NAMSET ( I ) ) GO TO 25 
20 CONTINUE 

LEVOLD = LEVEL 
LEVEL' = IER 
IER = LEVOLD 

IF (LEVEL. LT.O) LEVEL = 4 
IF (LEVEL. GT. 4) LEVEL = 4 
GO TO 30 
25 CONTINUE 

IF (LEVEL.LT. 4) GO TO 30 
CALL UGETI0(1,NIN, IOUNIT) 

IF (IEQDF.EQ. 1) WRITE (IOUNIT, 50) IER .NAMEQ, IEQ, NAME 
IF (IEQDF.EQ. 0) WRITE (IOUNIT, 50) IER.NAME 
30 IEQDF = 0 • 

RETURN 

35 FORMAT (19H *** TERMINAL ERROR, 10X, 7H( IER = ,13, 

1 20H) FROM IMSL ROUTINE ,3A2,A1,3A2) 

40 FORMAT ( 3 6H *** WARNING WITH FIX ERROR (IER = ,13, 

1 2 OH) FROM IMSL ROUTINE ,3A2,A1,3A2) 

45 FORMAT (18H *** WARNING ERROR, 1 IX, 7H( IER = ,13, 

1 20H) FROM IMSL ROUTINE ,3A2,A1,3A2) 

50 FORMAT (20H *** UNDEFINED ERROR, 9X, 7H( IER = ,15, 

1 20H) FROM IMSL ROUTINE ,3A2,A1,3A2) 

55 IEQDF = 1 
DO 60 1=1,3 
60 NAMEQ(I) = NAME (I) 

65 RETURN 
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ZSC03680 
ZSC03690 
ZSC03700 
ZSC03710 
ZSC03720 
ZSC03730 
ZSC03740 
ZSC03750 
ZSC03760 
ZSC03770 
ZSC03780 
ZSC03790 
ZSC03800 
ZSC03810 
ZSC03820 
ZSC03830 
ZSC03840 
ZSC03850 
ZSC03860 
ZSC03870 
ZSC03880 
ZSC03890 
ZSC03900 
ZSC039 10 
ZSC03920 
ZSC03930 
ZSC03940 
ZSC03950 
ZSC03960 
ZSC03970 
ZSC03980 
ZSC03990 
ZSC04000 
ZSC04010 
ZSC04020 
ZSC04030 
ZSC04040 
ZSC04050 
ZSC04060 
ZSC04070 
ZSC04080 
ZSC04090 
ZSC04100 
ZSC04110 
ZSC04120 
ZSC04130 
ZSC04140 
ZSC04150 
ZSC04160 
ZSC04170 
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END ZSC04180 

ZSC0419Q 
ZSC04200 
ZSC04210 
ZSC04220 
ZSC04230 
ZSC04240 
ZSC04250 
ZSC04260 
ZSC04270 
ZSC04280 
ZSC04290 
ZSC04300 
ZSC04310 
ZSC04320 
ZSC04330 


SUBROUTINE VXADD(A.ACC) 

DOUBLE PRECISION A, ACC (2) 

DOUBLE PRECISION X,Y,Z,ZZ 

X = ACC(l) 

Y = A 

IF (DABS (ACC ( 1 ) ) . GE . DABS (A) ) GO TO 1 
X = A 

Y = ACC(l) 

1 Z = X+Y 

ZZ = (X-Z)+Y 
ZZ = ZZ+ACC(2) 

ACC(l) = Z+ZZ 

ACC (2) = (Z-ACC(1))+ZZ 

RETURN 

END 


ZSC04340 

ZSC04350 

ZSC04360 

ZSC04370 

ZSC04380 

ZSC04390 

ZSC04400 

ZSC04410 

ZSC04420 

ZSC04430 

ZSC04440 

ZSC04450 

ZSC04460 

ZSC04470 

ZSC04480 


SUBROUTINE UGETIO(IOPT,NIN,NOUT) 


INTEGER 
INTEGER 
DATA 

IF (IOPT.EQ.3) GO TO 10 

IF (IOPT.EQ.2) GO TO 5 

IF (IOPT.NE.l) GO TO 9005 

NIN = NIND 
NOUT = NOUTD 
GO TO 9005 
5 NIND = NIN 
GO TO 9005 
10 NOUTD = NOUT 
9005 RETURN 
END 


IOPT, NIN, NOUT 
NIND, NOUTD 
NIND/5/.N0UTD/6/ 
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SUBROUTINE VXMUL l 

(A.B.ACC) 

ZSC04490 

DOUBLE PRECISION 

A,B , ACC(2) 

ZSC04500 

DOUBLE PRECISION 

X,HA,TA,HB,TB 

ZSC04510 

INTEGER 

IX(2) , I 

ZSC04520 

LOGICAL* 1 

LX(8) ,LI(4) 

ZSC04530 

EQUIVALENCE 

(X,LX(1),IX(1)) 4 (I,LI(D) 

ZSC04540 

DATA 

1/0/ 

ZSC04550 

X = A 


ZSC04560 

LI (4) = LX (5) 


ZSC04570 

IX ( 2 ) = 0 


ZSC04580 

I = (I/16)*16 


ZSC04590 

LX(5) = LI (4) 


ZSC04600 

HA=X 


ZSC04610 

TA=A-HA 


. ZSC04620 

X = B 


ZSC04630 

LI (4) = LX(5) 


ZSC04640 

IX(2) = 0 


ZSC04650 

I = (1/16 )* 16 


ZSC04660 

LX(5) = LI (4) 


ZSC04670 

HB = X 


ZSC04680 

TB = B-HB 


ZSC04690 

X = TA*TB 


ZSC04700 

CALL VXADD(X.ACC) 


ZSC04710 

X = HA*TB 


ZSC04720 

CALL VXADD(X, ACC) 


ZSC04730 

X = TA*HB 


ZSC04740 

CALL VXADD(X,ACC) 


ZSC04750 

X = HA*HB 


ZSC04760 

CALL VXADD(X,ACC) 


ZSC04770 

RETURN 


ZSC04780 

END 


ZSC04790 

SUBROUTINE VXSTO 

(ACC ,D) 

ZSC04800 

DOUBLE PRECISION 

ACC (2) ,D 

ZSC04810 

D = ACC(1)+ACC(2) 


ZSC04820 

RETURN 


ZSC04830 

END 


ZSC04840 


177 r 




FILE: ZSCNT FORTRAN A1 IMSL FRI 06/17/83 08:20:26 


PAGE 12 OF 14 


SUBROUTINE ZSCNU (X,N,FCN,NDIGIT,N1,A,Z,Y,XN0RM,B,WK, 
1 MAXIT.PAR, IER) 

INTEGER 

DOUBLE PRECISION 

1 

INTEGER 

1 

DOUBLE PRECISION 

1 

DOUBLE PRECISION 
DATA 
IER = 0 

DSEED = 12345. ODO 
CFACT = 0.99D0 
BIG = 5.0D5 
SFACT = 0 . IDO 
IBMAX = 50 
NRS = 2 

RACC = DMIN1 (DMAX1 (SMALL, 10 . ODO**(-NDIGIT) ) , 0 . IDO) 

REPS = DSQRT (SMALL) 

HLMAX = 3. ODO 
ITER = 0 
IEVAL = 0 
NSTART = 0 
IBNORM = 0 
RX = l.ODO 
RRX = O.ODO 
EPS = O.ODO 
DO 5 1=1, N 
B (10 = 0 . DO 
A(N1 , I) = l.ODO 
Z (I , N1 ) = X(I ) 

5 CONTINUE 
B(N1) = 1 • DO 
A(N1 ,N1) = l.ODO 
JI = N1 

IEVAL = IEVAL+1 

CALL FCN (Z ( 1 ,N1) ,A(1,N1) ,N, PAR) 

BNORM = O.ODO 
DO 10 1=1, N 

BNORM = BNORM+A ( I , N 1 ) *A ( I , N 1 ) 

10 CONTINUE 

XNORM(Nl) = BNORM 

15 IF (NSTART. EQ. NRS) EPS = EPS* 10. ODO 
IF (NSTART. EQ.O) EPS = DMIN1(RX,REPS) 

IF (EPS. GT. BIG) GO TO 120 
NSTART = MOD (NSTART, NRS )+l 
DO 30 J=1,N 
DO 25 1=1, N 

20 TR = (GGUBFS (DSEED) -0.5D0)*2. ODO 

IF (DABS (TR ) .LT.O. IDO) GO TO 20 
Z(I,J) = X(I)+DMAX1 (DABS (X(I ) ) , 0 . 1D0)*TR*EPS 
25 CONTINUE 

IEVAL = IEVAL+1 

CALL FCN (Z(1,J),A(1,J), N, PAR) 


IER,MAXIT,N,NDIGIT,N1 

A(N1,1),B(1), PAR( 1) ,WK( 1) ,X( 1) ,XN0RM(1) , 

Y(l) ,Z(N, 1) 

I, IBMAX, IBNORM,. IDGT, IEVAL, ITER, J,JER,JI,JS, 
MI, NRS, NSTART 

BIG, BNORM , CFACT , DX , EPS , HALF , HLMAX , RACC , REPS , 

RRX , RX , SFACT , SMALL , TEST , TN , TR 

DSEED 

SMALL/ Z34 100000 00000000/ 


ZSC04850 

ZSC04860 

ZSC04870 

ZSC04880 

ZSC04890 

ZSC04900 

ZSC049 10 

ZSC04920 

ZSC04930 

ZSC04940 

ZSC04950 

ZSC04960 

ZSC04970 

ZSC04980 

ZSC04990 

ZSC05000 

ZSC05010 

ZSC05020 

ZSC05030 

ZSC05040 

ZSC05050 

ZSC05060 

ZSC05070 

ZSC05080 

ZSC05090 

ZSC05100 

ZSC05110 

ZSC05120 

ZSC05130 

ZSC05140 

ZSC05150 

ZSC05160 

ZSC05170 

ZSC05180 

ZSC05190 

ZSC05200 

ZSC05210 

ZSC05220 

ZSC05230 

ZSC05240 

ZSC05250 

ZSC05260 

ZSC05270 

ZSC05280 

ZSC05290 

ZSC05300 

ZSC05310 

ZSC05320 

ZSC05330 

ZSC05340 

ZSC05350 

ZSC05360- 

ZSC05370 

ZSC05380 

ZSC05390 
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30 CONTINUE 
DO 35 J=1 , N 

XNORM(J) = O.DO 
DO 35 1=1, N 

XNORM(J) = XNORM(J)+A(I,J)*A(I,J) 

35 CONTINUE 
40 JI = N1 
JS = JI 
DO 45 J=1 , N 

IF (XNORM(J) .GT.XNORM(JS) ) JS = J 
IF (XNORM(J) . LT.XNORM(JI) ) JI = J 
45 CONTINUE 

IF (XNORM(JI) .EQ.O.DO) GO TO 125 
IF (XNORM(JI) .GT.SFACT*BNORM) GO TO 50 
BNORM = XNORM(JI) 

I BNORM = ITER 

50 IF ( (ITER -I BNORM) . GT. IBMAX) GO TO 120 
ITER = ITER+1 

IF (ITER.GE.MAXIT) GO TO 115 
DO 55 MI=1 ,N1 
. Y(MI) = B(MI) 

55 CONTINUE 
IDGT = 0 

CALL LEQT2F (A, 1,N1,N1,Y,IDGT,WK,JER) 

IF (JER.NE.O) GO TO 85 
DO 65 1=1, N 
DX = O.ODO 
DO 60 J=1,N1 

DX = DX + Y(J)*Z(I, J) 

60 CONTINUE 
X(I) = DX 
65 CONTINUE 
HALF = O.DO 
70 IEVAL = IEVAL+1 

CALL FCN (X,Y,N,PAR) 

TN = O.DO 
DO 75 1=1, N 

TN = TN+Y(I)*Y(I) 

75 CONTINUE 

IF (TN.LT.XNORM(JS) ) GO TO 95 
HALF = HALF+1 . DO 
IF (HALF . GT . HLMAX) GO TO 85 
DO 80 1=1, N • 

X(I) = (X(I)+HALF*Z (I , JI ))/ (HALF+1 . ODO) 

80 CONTINUE 
GO TO 70 

85 IF (JI.EQ.N1) GO TO 15 
XNORM(Nl) = XNORM(JI) 

DO 90 1=1, N 

Z(I,N1) = Z ( I , JI) 

A(I ,N1 ) = A ( I , JI) 

90 CONTINUE 
GO TO 15 

95 IF ((HALF. NE. O.DO). OR. (ITER. EQ.l)) GO TO 105 
RX = SMALL 


ZSC05400 

ZSC05410 

ZSC05420 

ZSC05430 

ZSC05440 

ZSC05450 

ZSC05460 

ZSC05470 

ZSC05480 

ZSC05490 

ZSC05500 

ZSC05510 

ZSC05520 

ZSC05530 

ZSC05540 

ZSC05550 

ZSC05560 

ZSC05570 

ZSC05580 

ZSC05590 

ZSC05600 

ZSC05610 

ZSC05620 

ZSC05630 

ZSC05640 

ZSC05650 

ZSC05660 

ZSC05670 

ZSC05680 

ZSC05690 

ZSC05700 

ZSC05710 

ZSC05720 

ZSC05730 

ZSC05740 

ZSC05750 

ZSC05760 

ZSC05770 

ZSC05780 

ZSC05790 

ZSC05800 

ZSC05810 

ZSC05820 

ZSC05830 

ZSC05840 

ZSC05850 

ZSC05860 

ZSC05870 

ZSC05880 

ZSC05890 

ZSC05900 

ZSC05910 

ZSC05920 

ZSC05930 

ZSC05940 
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DO 100 1=1, N 

ZSC05950 


RX = DMAX1 (RX,DABS(X(I ) -Z(I , JI) )/DMAXl (DABS (X(I) ) , 0 . IDO) ) 

ZSC05960 

100 

CONTINUE 

ZSC05970 


RRX = DMAXI(-DLOGIOCRX) ,0.0D0) 

ZSC05980 


IF (RX.LE.RACC) GO TO 125 

ZSC05990 

105 

IF (TN.LT.CFACT*XNORM(JI)) N START = 0 

ZSC06000 


XNORM(JS) = TN 

ZSC06010 


DO 110 1=1, N 

ZSC06020 


Z(I,JS) = X(I) 

ZSC06030 


A(I , JS) = Y C I ) 

ZSC06040 

110 

CONTINUE 

ZSC06050 


GO TO 40 

ZSC06060 

115 

IER = 129 

ZSC06070 


GO TO 125 

ZSC06080 

120 

IER = 130 

ZSC06090 

125 

DO 130 1=1, N 

ZSC06100 


X(I) = Z(I , JI) 

ZSC06110 

130 

CONTINUE 

ZSC06120 


RETURN 

ZSC06130 


END 

ZSC06140 
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APPENDIX C 


RESULTS OP EXPERIMENTS - NEW SERIES 



_J I i 1 I I I I I I I I I L 

95 100 105 110 115 120 125 130 135 140 145 150 155 160 

Temperature (°F) 

Fig. C-l Viscosity of 20W40 Oil 
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TABLE D-l 


PARAMETERS USED IN CALCULATIONS 


RUN 

MATERIAL 

R 

mm ( in . ) 

c-io 3 

mm ( in . ) 

t 

mm (in. ) 

E a ' 
MPa-10 

(psi- 10^) 

V 

L 1 

mm ( in. ) 

L 

mm (in.) 

e 

mm ( in . ) 

I 

Babbitt 

9.52 

(0.375) 

11.43 

(0.45) 

1.19 

(.047) 

51.7 

(7.5) 

0.36 

7.56 

(0.298) 

6.78 

(0.267) 

2.92 

(.115) 

II 
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9.52 

(0.375) 

19.05 

(0.45) 

1.19 

(.047) 
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0.36 

7.56 

(0.298) 

6.78 

(0.267) 

2.92 

(.115) 

V 

Rulon J 

9.52 

(0.375) 

42.55 
(1.67 5) 

2.44 

(.096) 

1.72 

(.25) 

0.46 

7.56 

(0.298) 

6.78 

(0.267) 

2.92 

(.115) 

VI 

Carbon 

9.52 

(0.375) 

20.35 

(0.80) 

1.44 

(.056) 

21.5 

(3.11) 

0.29 

7.56 

(0.298) 

6.78 

(0.267) 

2.92 

(.115) 

VII 

Babbitt 

6.00 

(0.236) 

8.89 

(0.35) 

0.89 

(.035) 

51.7 

(7.5) 

0.36 

5.59 

(0.22) 

5.08 

(0.2) 

1.93 

(.076) 


n 

o 

73 

73 

m 
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o 

H 

X 
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O 

po 

K 


*The viscosity used is that corresponding to test temperature plus an assumed 20°F used in the oil film. 


FIGURES SHOWING COMPARISON WITH EXPERIMENTAL DATA 
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Fig. D-4 Performance of Babbitt Pumping Ring (I) 
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Fig, D-6 Performance 
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Fig. D-7 Performance of Babbitt Pumping Ring (II) 
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Fig. D-9 Performance of Small Babbitt Pumping Ring (VII) 
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Fig. D-10 Performance of 


Stroke = 50.8 mm 
p 0 = 6.89 MPa 
C =8.90 micron 

o = 10 Hz Test 

A = 35 Hz Test 

x = 60 Hz Test 

= Theory 
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Fig. D-ll Performance of Small Babbitt Pumping Ring (VII) 
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Fig. d- 13 Qualitative Comparisons of Theoretical and Experimental Results 
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Fig. D-15 Performance of Carbon Pumping Ring 
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Fig. 0-16 Performance of Carbon Graphite Pumping Ring 
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Fig. D-20 Performance of Rulon Pumping Ring 
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TABLE D-2 


QUALITATIVE SUMMARY OF COMPARISON 
BETWEEN ANALYTICAL RESULTS AND TESTS 


RING 

THE0 /Q o EXP^ 

(p fm THE0 /p fm EXP ) 

H igh p Q 

Low p 
r o 

High p Q 

Low p 

o 

Babbitt 

Small C 

2 

1 

2 

1 

Large C 

3 

1.5 

BBS 

0.7 - 0.2 

Small R 

2 

0.5 

2 

IB 

Carbon Graphite 

2.5 

3-2 

1-2 

1.5 - 4 

Rulon 

2.5 

2.5 

1 

BBS 
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